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By asymptotically matching a post-Newtonian (PN) metric to two perturbed Schwarzschild met- 
rics, we generate approximate initial data (in the form of an approximate 4-metric) for a nonspinning 
black hole binary in a circular orbit. We carry out this matching through 0(v 4 ) in the binary's or- 
bital velocity v, and thus the resulting data, like the 0(v 4 ) PN metric, are conformally curved. 
The matching procedure also fixes the quadrupole and octupole tidal deformations of the holes, 
including the 1PN corrections to the quadrupole fields. Far from the holes, we use the appropriate 
PN metric that accounts for retardation, which we construct using the highest-order PN expressions 
available to compute the binary's past history. The data set's uncontrolled remainders are thus 
0(v 5 ) throughout the timeslice; we also generate an extension to the data set that has uncontrolled 
remainders of 0{v 6 ) in the purely PN portion of the timeslice (i.e., not too close to the holes). 
This extension also includes various other readily available higher-order terms. The addition of 
these terms decreases the constraint violations in certain regions, even though it does not increase 
the data's formal accuracy. The resulting data are smooth, since we join all the metrics together 
by smoothly interpolating between them. We perform this interpolation using transition functions 
constructed to avoid introducing excessive additional constraint violations. Due to their inclusion of 
tidal deformations and outgoing radiation, these data should substantially reduce both the high- and 
low-frequency components of the initial spurious ( "junk" ) radiation observed in current simulations 
that use conformally flat initial data. Such reductions in the nonphysical components of the initial 
data will be necessary for simulations to achieve the accuracy required to supply Advanced LIGO 
and LISA with the templates necessary for parameter estimation. 

PACS numbers: 04.25.Nx, 04.25. dg, 04.30.Db 



I. INTRODUCTION 
A. Motivation and overview of results 

At present, several years after the initial breakthroughs 
in the evolution of binary black hole spacetimes [l], H, Q , 
numerical relativity has matured to the point where suc- 
cessful binary black hole simulations are now common- 
place (see, e.g., Q for a review). Recent progress in the 
nonspinning case includes simulations of systems with 
mass ratios of up to 10 : 1 0, 0] and longer, more accu- 
rate simulations of equal-mass systems Q H| . 

It is now time to consider what improvements need to 
be made to these simulations so that they are accurate 
enough to provide gravitational wave detectors such as 
LIGO 0] with the model waveforms they need to detect 
and study binary black holes. The accuracy required of 
such waveforms has been studied by Lindblom, Owen, 
and Brown (LOB) [l0|. Their results imply that cur- 
rent simulations are sufficiently accurate to supply the 
waveforms necessary for detection with either LIGO or 
LIS A uA\ . This conclusion is supported by the Samu- 
rai [l2| and NINJA [l3| projects, which indicate that 
currently used data analysis pipelines (including some 
not based on matched filtering) can easily detect a wide 
variety of numerical relativity waveforms at essentially 



the same level in stationary gaussian noise. The Samurai 
project also performs a more detailed comparison of a 
subset of waveforms for parameter estimation, and finds 
that they are all indistinguishable if used for estimation 
of intrinsic parameters (i.e., not sky position or arrival 
time) at a signal-to-noise ratio (SNR) of less than 14 (25 
if one eliminates the code that disagrees the most with 
the others). 

However, Advanced LIGO may detect binary black 
hole signals with SNRs of order 100 [l(|. According 
to LOB [n|, if a waveform's phase error (suitably av- 
eraged in the frequency domain) is less than 0.007 radi- 
ans, then it is suitable for use in parameter estimation 
with Advanced LIGO with such an SNR. However, even 
the Caltech/Cornell group's simulations (which are ar- 
guably the most accurate yet) have maximum phase er- 
rors (in the time domain) of order 0.01 radians or more 
(see 0, for some discussion of their error budget). 
Converting such error measures to the criteria of LOB 
is more subtle than it looks and we do not attempt it 
here. (See [l5[ for a discussion of some possible pitfalls, 
along with suggestions for successful applications of the 
standards from LOB.) However it seems likely that even 
the Caltech/Cornell group's simulations may not satisfy 
the LOB conditions for parameter estimation, at least for 
binaries whose masses place the worst phase error in the 
detector's most sensitive frequency band. 
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All current binary black hole simulations incur some 
error from the initial data used. These data sets' lack of 
astrophysical realism is clearly announced by the burst of 
spurious (or "junk") radiation present at the beginning 
of these simulations. This junk radiation is also respon- 
sible for various deleterious effects on simulations. First, 
one wastes the computer time required for the spurious 
radiation to propagate off the computational grid. Sec- 
ond, the system that remains after that time has been 
perturbed by the spurious radiation. This radiation's 
most prominent effect is to increase the binary's eccen- 
tricity, 1 though it also slightly increases the masses of 
the holes — see [l9j . Additionally, in the unequal-mass 
case, the initial junk radiation is emitted anisotropically, 
giving the system a small "kick" transverse to the holes' 
initial orbital motion [2p| . The high-frequency compo- 
nent of the spurious radiation is a particular problem for 
spectral codes. For instance, the Caltech/Cornell group 
finds that the high-frequency component of the initial 
pulse of spurious radiation generates secondary spurious 
waves that propagate throughout the computational do- 
main for two light-crossing times after the initial junk 
radiation has exited @, [H| • 

Currently employed initial data's omission of signifi- 
cant features of the spacetime can also be seen analyti- 
cally. Except for occasional tests of initial data (discussed 
in Sec. II Bp . all current simulations use initial data sets 
that assume conformal flatness (i.e., that the spatial met- 
ric is a scalar multiple of the flat 3-metric). The assump- 
tion of conformal flatness is convenient, since it allows 
one to get simple, mostly analytic expressions for initial 
data that exactly solve the constraint equations and in- 
clude orbiting black holes (see, e.g., [2l|, [22j for reviews). 

In general, these sets are geared either towards the 
puncture or excision approaches. The majority of the 
community uses punctures, with initial data stemming 
from [23| . These data are very flexible, as they contain 
parameters with which one can directly set the momen- 
tum and spin of each hole. For instance, for evolutions 
with spinning holes, one can simply set these parameters 
using post-Newtonian (PN) results as in [2I I2H I2H . 
For nonspinning configurations in a circular orbit, pa- 
rameter choices based on the assumption of a helical 
Killing vector [28|, [2^] are possible as well. The Cal- 
tech/Cornell and Princeton groups use excision (30l. [3ll]. 
with initial data constructed using the conformal thin- 
sandwich method (see, e.g., [13]). These data are slightly 
harder to construct than puncture data are, since one has 
to solve a larger number of elliptic equations. However, 
excision data have the advantage of a direct connection to 
the isolated horizon formalism [33|, which allows one to 
construct holes with well-defined masses and spins. Ad- 
ditionally, the excision approach is applicable to a wider 



1 In practice, the eccentricity can be reduced by various means [l6t 
mill- 



array of initial data construction methods: It is used in 
all of the extant evolutions of superposed black hole data 
sets except for one specifically tailored to the puncture 
approach. (See Sec. II Bl for a discussion of these evolu- 
tions.) The data we present here also require evolutions 
using excision or the turducken approach 34 1. 

Conformally flat initial data cannot accurately repre- 
sent some features of a binary black hole spacetime, since 
the PN metric for a binary system stops being spatially 
conformally flat at 0(u 4 ), where v is the binary's orbital 
velocity in units of c, the speed of light (see Sec. IIV Al 
for discussion). This is the same order at which gravita- 
tional radiation enters the PN metric (see Sec. IIV Bl for 
discussion). The order at which this fundamental dis- 
agreement with PN predictions first occurs gives a rough 
indication of the error committed in using conformally 
flat initial data. At present, the simulation for which this 
error is the smallest is the longest of the Caltech/Cornell 
runs (in Q), for which initial - 0.24, where iWiai is 
the binary's initial orbital velocity. We expect the ini- 
tial data's conformal flatness to only affect the waveform 
at 0(vf nitial ), which for this run is comparable to the 
phase error allowed for a waveform to be used for param- 
eter estimation with Advanced LIGO. It is thus possible 
(though perhaps not likely) that conformally flat initial 
data would be suitable for use in the simulations that 
will generate such waveforms. 

It is unlikely that conformally flat initial data can be 
used to generate the waveforms required for parameter 
estimation with LISA. Here the required (appropriately 
averaged frequency domain) phase accuracy is 2 x 
radians [l(| , 20 times smaller than vf nitiai for the longest 
of Caltech/Cornell's simulations. One can reduce the er- 
ror in the initial data by starting the simulation with 
a larger separation. However, f-^itiai — 2 x 10~ 4 im- 
plies an initial (PN coordinate) separation of ~ 71 times 
the binary's mass, and thus a merger time that is over 
400 times as long as the Caltech/Cornell group's longest 
simulations to date, which start from a (PN coordinate) 
separation of ~ 15.3 times the binary's mass. It is thus 
necessary to improve the accuracy of the initial data. 
Evolutions of more accurate initial data will also give a 
direct measure of the errors introduced in using current, 
conformally flat initial data. 

This paper provides initial data that include more of 
the physics present in the binary's spacetime than any 
previous constructions. In particular, our data's accu- 
rate description of certain properties of the spacetime 
should substantially reduce both components of the afore- 
mentioned spurious radiation. These two components 
are thought to come from different physical effects. The 
long-wavelength component is thought to correspond to 
the initial data's lack of outgoing gravitational radiation, 
whose wavelength would be somewhat longer than the 
orbital separation. Of course, one expects the junk to be 
generated predominantly in the strong-field region near 
the holes, where the binary's gravitational radiation can- 
not be disentangled from the rest of its gravitational field. 
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However, one also expects the pieces that one wants in 
the strong-field region to appear at 0(v 4 ), just as the true 
gravitational waves do in the radiation (or far) zone (de- 
fined in Sec. [TT|- Our data include all the 0(v 4 ) terms in 
the strong and weak field regions. The short- wavelength 
component is thought to come from the holes' quasinor- 
mal modes ringing down, emitting gravitational radiation 
with wavelengths on the order of their masses, as they 
relax from their initial, close to spherical state to their 
desired tidally deformed state (see, e.g., fl4{|). Our data 
include the Newtonian quadrupole and octupolc tidal de- 
formations each hole induces on the other, as well as the 
1PN corrections to the quadrupole deformations. 

The tidal deformations are contained in perturbed 
Schwarzschild metrics, given (in the horizon-penetrating 
coordinate system we use) in Eqs. fl3.5[) and (|3.6p . The 
tidal fields are fixed by asymptotically matching these 
Schwarzschild metrics to an 0(v 4 ) PN metric, given in 
Eqs. (|4.1|) ; expressions for the tidal fields around hole 1 
are given in Eqs. (|Bip . One also needs to introduce a 
coordinate transformation in order to put the black hole 
metrics in the same coordinate system as the PN met- 
ric. This transformation is also determined (perturba- 
tively) by the matching; instructions for putting it to- 
gether around hole 1 are given in Sec. IV Gl Instructions 
for converting all of these results to the region around 
hole 2 are given at the beginning of Sec. 

The PN metric mentioned above treats retardation 
perturbatively and thus becomes inaccurate far from the 
binary. In that region, we thus use a version of the PN 
metric that includes retardation explicitly (but also uses 
a multipolar decomposition, so it does not provide the de- 
sired accuracy closer to the holes). This metric is given 
in Eqs. (|6.6p . Due to retardation, one needs to know the 
binary's past history accurately in order to obtain the 
far zone metric accurately far from the binary. This past 
history is computed in Sec. IVI Al to the highest PN or- 
der possible with current results. The contributions of 
these terms are of equal or higher order than some of our 
uncontrolled remainders, but their inclusion is necessary 
if one wishes to obtain, e.g., the correct phasing for the 
outgoing radiation. 

We have also added other formally higher-order terms 
to the metrics, including all the readily available PN re- 
sults, along with a resummation of the black hole back- 
grounds in the PN metric. We found that certain of 
these terms improved the constraint violations in various 
regions; other of these additional terms arc expected to 
improve evolutions of the data. The resummation is given 
in Sec. IVIII Al all the remaining higher-order terms are 
discussed in Sec. IVIII The transition functions that we 
use to stitch the metrics together smoothly are given in 
Sec. IVIII B I We have constructed these transition func- 
tions so that they satisfy the so-called Frankenstein the- 
orems (35j . The resulting merged metric is thus guaran- 
teed to have constraint violations whose formal order is 
no larger than the constraint violations of the individ- 
ual metrics. See Appendix [D] for the technical details of 
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TABLE I: The L 2 and sup norms (denoted by 1 1 ■ 1 1 2 and 1 1 ■ 1 1 oc , 
respectively) of the Hamiltonian and momentum constraints 
(H and M) for an equal-mass binary with a coordinate separa- 
tion of 10 times its total mass, m. (The norms of the momen- 
tum constraint include the L 2 3- vector norm \/M k Mk, where 
the index is lowered using the metric under consideration.) 
These are computed along the x-axis outside the unperturbed 
horizons of the holes, but inside the interval [—16.4m, 16.4m]. 



how we compute the metrics. We have generated Maple 
scripts and C code that produce initial data for a non- 
spinning binary of any mass ratio and initial separation 
in a quasicircular orbit. These are available at [361 ] - 
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FIG. 1: The Hamiltonian constraint and norm of the momen- 
tum constraint along the ic-axis around hole 1 for this pa- 
per's data along with the data from Papers I and II. (In the 
norm of the momentum constraint, the index is lowered using 
the metric in question.) All of these were computed for an 
equal-mass binary with a coordinate separation of 10 times its 
mass. In the inset, we zoom in to show how the Hamiltonian 
constraint violations behave close to the horizon. Note that 
the data from Paper II are in a different coordinate system 
throughout, and that the black hole background in Paper I's 
data is in a different (and not horizon penetrating) coordinate 
system. 

The constraint violations give a measure of how much 
more accurate this paper's data are c omp ared with those 
constructed in previous papers (viz., [37| and [38| . which 
we refer to as Papers I and II). We plot the constraint 
violations of the three sets of data in Fig. [TJ [See, e.g., 
Eqs. (14)-(15) in [2l[ for expressions for the constraint 
equations.] N.B.: While the Hamiltonian constraint vi- 
olations of this paper's data are larger close to the hole 
than those of cither of the previous papers' data, this is 
due to our inclusion of the full time dependence of the 
tidal fields, as discussed in Appendix ID II the motivation 
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for doing this, even at the cost of larger constraint viola- 
tions, is given in Sec. IVII1 If one does not include these 
higher-order terms, then the new data's constraint vio- 
lations are smaller close to the hole than those of either 
of the previous papers' data. See Fig. [20] for a compar- 
ison of the constraint violations of the inner zone met- 
ric with and without full time dependence in the tidal 
fields. For a simple, quantitative comparison of the over- 
all constraint violations, we can consider their L 2 and 
sup norms. These are presented for all three sets of data 
in Table U and show, as expected, that the new data 
have smaller overall constraint violations than either of 
the previous papers' data. Of course, the L 2 norm of 
the new data's Hamiltonian constraint violations is only 
very slightly less than that of the data from Paper II, 
but this is probably to be expected: See Sec. IVIII CI for 
some discussion of this, as well as further details about 
the comparison plot and table. 



B. Conformally curved initial data 

The problems with conformally flat data have inspired 
a variety of constructions of conformally curved binary 
black hole data, for which there are two general ap- 
proaches. One approach — the one we have taken — is pri- 
marily motivated by a desire for astrophysical realism, 
seeing the spurious radiation as an indicator of the fail- 
ure of conformally flat initial data to model the desired 
spacetime accurately enough. These constructions so far 
have restricted their attention to nonspinning binaries 
and used the PN approximation to include the binary's 
physics. (Of course, true astrophysical binaries are ex- 
pected to have significant spins, so the consideration of 
nonspinning binaries is merely a technical convenience 
appropriate for initial attempts at constructing astro- 
physically realistic data.) The other approach is primar- 
ily concerned with reducing the junk radiation in prac- 
tice with relatively simple choices for initial data (viz., 
a superposition of boosted black holes). This approach 
is often geared primarily towards spinning black holes, 
since the amount of spurious radiation increases with the 
spin of the holes. (This is due to the nonexistence of 
a conformally flat slicing even of an isolated Kerr black 
hole with nonzero spin — see [H, E(|.) 

Besides our present work, other constructions in the 
first category are this work's antecedents ([13, [H, [4l| . 
discussed in the next subsection), along with the ap- 
proaches of Nissanke [12] and Kelly et al. |43| • (The latter 
two only use the PN approximation in their construc- 
tion, so the resulting data cannot accurately describe the 
spacetime near the holes.) Nissanke, building on the work 
of Blanchet [3] (who constructed initial data for a head- 
on collision of initially stationary holes) , obtains explicit 
analytic expressions for the 3-metric and extrinsic curva- 
ture from the 2PN metric. (N.B.: These constructions 
only use the version of the PN metric that treats retar- 
dation pcrturbativcly, so the resulting data rapidly lose 



accuracy away from the binary.) Kelly et al. [43[ extend 
the work of Tichy et al. [45[ to give initial data from 
the 2.5PN ADMTT metric that are valid through 0(v 4 ) 
and thus contain the binary's outgoing radiation in the 
far zone. (See Sec. HVBI for a discussion of why gravita- 
tional radiation is present at that order.) They do not ob- 
tain data that are valid through 0(v 5 ), even though the 
2.5PN metric contains the 0(v 5 ) terms in all its com- 
ponents, since one would need the 0(v 6 ) terms in the 
spatiotemporal components to obtain initial data valid 
through 0(v 5 ): See Sec. [I C II for further details. In order 
to obtain these data, they have to evaluate an integral 
numerically, so their data are not completely analytic. 
The original data in Tichy et al. [45[ are completely ana- 
lytic, though they use the version of the ADMTT metric 
available in the literature, which gives expressions for the 
0{v A ) and 0(v 5 ) pieces of the transverse-traceless (TT) 
part of the metric that are only valid in the near zone. 
Tichy et a/.'s data were thus only valid through 0(v 3 ) 
far from the holes. Kelly et al. calculated the additional 
terms necessary for the TT pieces to be valid through 
0{v A ) in the far zone. 

The most recent progress in the second category for 
nonspinning black holes is Lovelace's [46| construction 
and evolution of superposed Schwarzschild data. (There 
is a companion construction and evolution of superposed 
spinning black hole data in Lovelace et al. [47j .) Lovelace 
uses the superposed black hole data as free data for a con- 
straint solver and finds that the resulting data produce 
less spurious radiation than conformally flat data. There 
are also older constructions in a similar vein [H, [49|, H3| , 
though their data have only been evolved in head-on col- 
lisions without first solving the constraints 51|. In this 
case, it was found that there was more spurious radiation 
with the superposed data than with Brill-Lindquist data. 
The Kerr puncture data presented by Hannam et al. [52j 
(who use Dain's work in [531 . [|54| to construct superposed 
Kerr free data in puncture form for a constraint solver) 
have been shown to reduce the junk radiation for head- 
on collisions of spinning black holes. However, the given 
construction is only applicable to the case of a head-on 
collision of initially stationary holes: It has not been ex- 
tended to give the holes linear momenta. The freedom 
inherent in these superposed black hole initial data con- 
structions has been studied in (55j . and a helical Killing 
vector version of superposed data is constructed in [56| . 
But neither of those papers' data sets have been evolved, 
to our knowledge. 

While the use of a conformally curved Kerr metric 
near each of the holes is likely responsible for most of 
the reduction of spurious radiation seen in the case of a 
spinning binary, it is not clear, physically, what feature 
of Lovelace's data is reducing the spurious radiation in 
the Schwarzschild case. Furthermore, it is known [42j 
that the superposed nonspinning Kerr-Schild metric dif- 
fers from the PN metric at 2PN, due to the lack of in- 
teraction terms (e.g., terms that involve products of the 
distances from the field point to each of the holes) . Such 
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a disagreement with the PN metric is probably charac- 
teristic of all superposed initial data constructions, thus 
limiting their ability to substantially reduce the junk ra- 
diation. 

Evolutions of conformally curved binary initial data 
are uncommon, and most have restricted their atten- 
tion to the computationally cheap case of head-on col- 
lisions, as seen above. Only Lovelace [4(| and Lovelace 
et al. have evolved superposed black hole data for 
an orbiting binary. For the PN data, such head-on tests 
are mostly inapplicable: Only Nissanke's data set is not 
already specialized to a (quasi)circular orbit. But even 
though her data (which reduce to Blanchet's in the limit 
of a head-on collision of initially stationary black holes) 
could be tested using a head-on collision, they have not 
been evolved, to our knowledge. 

In fact, the only evolution of PN initial data of which 
we are aware is that of Kelly et al. H3- They evolved 
the data they obtained in 1431] (as well as the original ver- 
sion without waves from |45j). necessarily doing so for an 
orbiting binary, since their data were derived under the 
assumption of a quasicircular orbit. They found what 
naively appears to be a slightly better reduction of the 
initial spurious radiation than that seen with Lovelace's 
superposed black hole initial data, even without first solv- 
ing the constraint equations. (Their data satisfy the con- 
straints approximately, but not exactly.) However, the 
results are not directly comparable: Most importantly, 
Lovelace begins his evolution with a separation of about 
twice that with which Kelly et al. begin their evolution. 
Thus, while Kelly et al. see a slightly larger reduction 
in the maximum amplitude of their junk radiation, even 
their smallest amplitude is larger than that from the con- 
formally flat data with which Lovelace is comparing his 
data's performance. In addition, Kelly et al. are compar- 
ing their data's performance with puncture data, while 
Lovelace is using conformal thin-sandwich excision data 
(constructed in the manner of Cook and Pfeiffer [S^]) as a 
benchmark. (There are also other important differences, 
such as the extraction radius used and the mode of the 
junk radiation that is reduced most substantially.) 

It is also important to realize that Kelly et al.'s data 
only reduce the low-frequency component of the junk ra- 
diation: The high-frequency component visible when us- 
ing conformally flat puncture data is still present in the 
evolution of Kelly et aZ.'s data. This is not surprising, 
since their initial data make no attempt to include accu- 
rate tidal deformations on the holes. (And, indeed, their 
data's constraint violations are largest near the holes.) 
In addition, they obtain nearly identical junk radiation 
when evolving the data with and without waves: The 
only difference is that the spurious radiation is super- 
posed over the outgoing wave train when they evolve the 
data with waves; the spurious radiation itself appears 
to be unchanged. This, again, is what should be antici- 
pated, since one expects the junk to be generated primar- 
ily in the strong field region near the holes, and Kelly et 
aL's data (with waves in the far zone) have the same ac- 



curacy in the strong field region as Tichy et a/.'s original 
data (without waves in the far zone). 

As mentioned in the previous subsection, the data set 
we construct here should help reduce both components 
of the junk radiation, since it includes an accurate de- 
scription of the spacetime near the holes (including the 
quadrupolar and octupolar tidal deformations) , matched 
to the PN metric through 0(v 4 ). We also offer an ex- 
tension to these data that are accurate through 0(v 5 ) in 
the PN portions of the timeslice (i.e., not too close to the 
holes) . This extension includes certain other higher-order 
terms, as well, including higher-order corrections to the 
trajectories. While these "extra" higher-order terms do 
not increase the data's formal accuracy, even in just the 
PN portions, such terms will likely improve the accuracy 
of the data in practice. This extension will also allow 
for a more direct comparison with Kelly et a/.'s results, 
since they included the 0(v 5 ) terms in the spatial metric 
in the near zone (though not the matching 0(v e ) terms 
in the extrinsic curvature, which we have). Such a com- 
parison would reveal how much of the spurious radiation 
is due to the initial data's failure to include the correct 
tidal deformations. 



C. Specifics of our approach and its relation to 
other work 

With currently available technology, if one wishes to 
generate initial data that include the holes' tidal defor- 
mations or the binary's outgoing radiation, it is neces- 
sary to allow the initial data (as first constructed) to be 
merely an approximate solution of the constraint equa- 
tions. (One can always use these approximate data as 
free data for a constraint solver and thus obtain an ex- 
act solution of the constraints, to numerical precision.) 
Since the post-Newtonian approximation has been devel- 
oped to a very high order, it is an obvious choice for the 
description of the binary's spacetime. Indeed, an explicit 
expression for the metric to 2.5PN order (with a per- 
turbative treatment of retardation) is given by Blanchct, 
Faye, and Ponsot (BFP) 

However, one cannot obtain accurate initial data 
throughout a timeslice of the binary's spacetime using 
just the PN metric, since the PN approximation breaks 
down near the holes: The PN approximation is a weak- 
field approximation (due to the post-Minkowskian itera- 
tion in powers of G, Newton's gravitational constant), in 
addition to being a slow-motion approximation (the post- 
Newtonian expansion proper, which formally proceeds as 
an expansion in 1/c, where c is the speed of light). More- 
over, the standard PN approximation (as presented, e.g., 
in BFP) treats retardation perturbatively. It thus be- 
comes inaccurate quite rapidly as one enters the radia- 
tion zone (i.e., when one is further than about a reduced 
gravitational wavelength away from the binary's center 
of mass). 

The resolution of both of these problems is to realize 
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that there is an appropriate approximate description of 
the spacetime in each of the regions where the standard 
PN metric is no longer a good approximation: Near each 
of the holes, in the regions known as the inner zones, 
spacetime is well described by a perturbed black hole 
metric. (These zones, as well as the others we mention 
here, are defined more precisely in Sec.|TT]) In the radia- 
tion zone (or far zone) there is another version of the PN 
metric that incorporates retardation nonperturbativcly. 
These are all readily available in the literature to the 
order we need them. For the reasons described below, 
we choose to use Detweiler's perturbed black hole met- 
ric in addition to BFP's PN metric. We also choose 
to put together the far zone metric following the recipe 
and ingredients supplied by Pati and Will [§3, EiH ■ (In 
addition, we use Blanchet's results [62j for the evolution 
of the binary's phase with radiation reaction: These sec- 
ular effects are important in obtaining the far zone metric 
accurately, due to its dependence upon retarded time — 
see Sec. EDS) 

One then has to stitch all these spacetimes — the far 
zone, inner zones, and the near zone where the standard 
PN metric is valid — together into one global approximate 
metric. (Of course, this metric will be global in space, but 
not in time — i.e., it will only be accurate in a temporal 
neighborhood of a timeslice of the binary's spacetime.) 
This stitching-together proceeds in two steps: First, one 
uses the technique of matched asymptotic expansions to 
match the metrics at a formal level. This puts all the met- 
rics in the same coordinate system (up to uncontrolled 
remainders) and fixes any previously undetermined pa- 
rameters (e.g., the holes' tidal perturbations) so that the 
metrics are asymptotic to each other in their regions of 
mutual validity (the buffer zones). The final, numeri- 
cal merging of the metrics is effected by transition func- 
tions that smoothly interpolate between the metrics in 
their mutual buffer zones. The resulting merged met- 
ric is guaranteed to satisfy the Einstein equations to the 
same order as its constituent metrics because the tran- 
sition functions are constructed following the so-called 
Frankenstein theorems [H[ . (We have checked this scal- 
ing explicitly in Sec. I VIII CI ) 

Once one has obtained such approximate initial data, 
it is, of course, possible to use them as the input to a 
numerical constraint solver, and thus obtain an exact so- 
lution, to numerical precision. In fact, it is probably 
desirable to do so. The idea when doing this is that if 
the input to the constraint solver satisfies the constraints 
to some reasonably good tolerance, and describes the 
desired physics, the "exact" solution one obtains after 
solving the constraints will not differ too much from the 
input in its physical content. This is probably true re- 
gardless of how one chooses to produce the "exact" solu- 
tion to the constraints, as long as the procedure modifies 
the initial guess in a reasonable way. (For instance, the 
York-Lichnerowicz decomposition [2l| multiplies the ini- 
tial data by an appropriate conformal factor and modifies 
the extrinsic curvature by the gradient of a vector field 



so that the data satisfy the constraint equations — both 
Pfeiffer et al. [63[ and Tichy et al. [4{| have implemented 
this numerically without the assumption of conformal 
flatness.) However, it may be preferable to project the 
approximate initial data onto the closest "point" on the 
constraint hypersurface (as measured by some appropri- 
ate norm) , possibly using the results of 64 1. 

If one chooses to evolve without solving the constraints 
(as did Kelly et al.), then one will not have a true vacuum 
evolution: The constraint violations will act as matter 
(which may not satisfy any of the standard energy con- 
ditions). For instance, Bode et al. [6{| investigated the 
evolution of initial data that only approximately satisfied 
the Hamiltonian constraint (though the momentum con- 
straint was satisfied exactly). They found that the holes 
accreted the negative Hamiltonian constraint violations 
that surrounded them, decreasing their masses. Addi- 
tionally, the initial apparent horizon masses were larger 
than the irreducible masses for the constraint-violating 
data, but equaled them for the constraint-satisfying data. 



Comparisons with similar constructions 



Alvi [41( was the first to attempt to construct bi- 
nary black hole initial data by matching perturbed black 
hole metrics onto a PN metric, considering, for sim- 
plicity, two nonspinning black holes in a circular orbit. 
(All subsequent attempts using this procedure, includ- 
ing the present one, have also restricted their attention 
to this simplest case.) In performing this calculation, 
he fixed the perturbation (encoded in multipolar tidal 
fields) on the holes a priori, using its expected Newto- 
nian quadrupole value, instead of reading it off from the 
matching as we do here. This gives the same result for 
the lowest-order pieces of the tidal fields as our method, 
though it docs not offer the opportunity of reading off 
the higher-order corrections, as we can. 

Numerical experiments with Alvi's initial data [66[ 
demonstrated that they were not suitable for use in evo- 
lutions. This inspired Yuncs et al. (37j (Paper I) to 
revisit the problem and correct various deficiencies in 
Alvi's method, such as an inconsistent order counting 
(due to not including the terms needed to compute the 
extrinsic curvature to the appropriate order; this is dis- 
cussed below) and a lack of smoothness at the joins of 
the global metric (due to not actually performing asymp- 
totic matching). Yunes and Tichy [38| (Paper II) then 
improved this initial data construction (in the sense of 
getting better numerical agreement between the metrics) 
by putting the near zone metric in a form that is very 
close to the inner zone metric near the holes. This was 
done by applying some resummation and using ADMTT 
coordinates instead of harmonic coordinates for the PN 
metric. They also constructed horizon-penetrating coor- 
dinates to give the first usable initial data generated with 
this method. However, in both cases the initial data were 
only valid through 0(v 2 ) and the tidal fields were still 
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just the lowest-order ones Alvi had obtained. 

The present calculation builds on all these previous 
attempts, computing fully matched initial data through 
0(v 4 ), so as to include the pieces of the PN metric that 
break conformal flatness (and contain gravitational radi- 
ation), and reading off the tidal fields (including the 1PN 
corrections to the quadrupole fields) from the match- 
ing. (We also demonstrate that the mass parameters 
of the PN metric are equal to those of the perturbed 
Schwarzschild metrics to the order considered: This is 
established here more firmly than in Papers I and II in 
addition to being extended to higher order.) In addition, 
we have included the radiation zone portion of the metric 
to accommodate the larger grid sizes common in current 
simulations. This was done explicitly by Alvi and implic- 
itly in Paper II, due to its use of the ADMTT PN metric, 
though neither of them included the effects of radiation 
reaction on the binary's past evolution in their far zone 
metric. The radiation zone was neglected completely in 
Paper I, which used a harmonic metric, as we do here. 
Additionally, all the previous versions used a corotating 
coordinate system, while the current calculation stays in 
inertial coordinates. 

While we only obtain fully matched initial data 
through 0(v 4 ), we actually have to carry out the match- 
ing of the 4- metric through 0(v 5 ) in order to do so: We 
need to match the 0(v 5 ) pieces of the spatiotcmporal 
components of the 4-metrics in order to obtain the extrin- 
sic curvature consistently [see the discussion after Eq. (2) 
in Paper I], and one needs to carry out the matching of 
all the components in order to obtain the 0(v 5 ) piece of 
the coordinate transformation. 

Additionally, while our goal was simply to keep terms 
of quadrupolar order overall in the multipole expansion, 
as was done previously, we found that it was necessary 
to match the lowest-order octupolar pieces in order to 
match the IPN corrections to the quadrupole pieces con- 
sistently. This is discussed in Sec. IV Al In fact, we 
have carried out the matching of quadrupole pieces to 
the highest possible order to which it can be done con- 
sistently without the inclusion of the hexadecapole tidal 
fields. (These hexadecapole pieces can only be included 
with input from nonlinear black hole perturbation theory, 
as the "quadrupole squared" pieces are of hexadecapole 
order.) We also obtained the IPN correction to the elec- 
tric octupole and the associated piece of the coordinate 
transformation as a further application of our matching 
procedure. However, we cannot obtain the other 0(v 4 ) 
octupole pieces in the initial data (since they include the 
hexadecapole tidal fields), so our knowledge of these cor- 
rections does not allow us to increase the formal order to 
which our data are valid. 

Building on the work done in Paper II, we have used 
horizon-penetrating coordinates for the black holes from 
the outset. This requirement of horizon penetration is 
necessary for numerical purposes. The coordinates need 
to be regular and the lapse positive in a neighborhood of 
the horizon: Even though the spacetime near the singu- 



larity will be excised or filled with matter, one needs to 
be able to evolve at least a small portion of the spacetime 
inside the horizon. 

At the same time, we want the coordinates for the 
black hole and PN metric to agree as closely as possible 
before the matching has been performed: Close agree- 
ment makes for simple matching algebraically and im- 
proves the numerical agreement of the resulting matched 
metrics. Ideally, the coordinates would agree exactly for 
an unperturbed black hole, though this is not compatible 
with the requirement of horizon penetration, as standard 
PN coordinates (harmonic or ADMTT) are not horizon 
penetrating. We thus attempted instead to obtain agree- 
ment between the two coordinate systems to as high a 
PN order as possible. 

These desiderata are satisfied if we use the fully har- 
monic version of Cook-Scheel coordinates (6?| for the 
black hole and standard (PN) harmonic coordinates for 
the PN metric: 2 Cook-Scheel coordinates are horizon- 
penetrating, and in their fully harmonic version only dif- 
fer from PN harmonic coordinates for an unperturbed 
black hole at 0(t> 4 ). Sec Appendix [A] for an explicit 
comparison. This agreement was the best of any of the 
horizon-penetrating coordinate systems present in the lit- 
erature we consulted (even if we also consider ADMTT 
coordinates for the PN metric). Of course, we then adjust 
this coordinate system perturbativcly so that it agrees 
with the near zone coordinate system to the order we 
have matched. However, we have checked that this ad- 
justment does not affect the coordinates' horizon pene- 
tration. The agreement between the coordinate systems 
used in Paper II was exact in the unperturbed case, while 
it is not here. However, we decided against converting 
the PN metric to horizon-penetrating coordinates for this 
version of the data — see Sec. IVIlfAl for further details. 

The other choices for our ingredients were made for 
computational ease. We selected Detweiler's perturbed 
black hole metric [5^] instead of Poisson's [H[ because 
Detweiler expresses the tidal fields in the Thorne-Hartle- 
Zhang (THZ) harmonic specialization of locally inertial 
coordinates [H, [7(J. This gauge choice agrees better 
with the PN metric in harmonic coordinates than does 
Poisson's light-cone gauge. For the far zone, the results 
from the direct integration of the relaxed Einstein equa- 
tions (DIRE) approach were the obvious choice: Pati and 
Will [10, Hl| give an explicit recipe for computing the far 
zone metric to the order we need it, along with all the nec- 



2 We use the term harmonic coordinates to refer to any coor- 
dinates x a that satisfy V a V a x' 3 = 0, not just PN harmonic 
coordinates. (Here V Q is the covariant derivative associated 
with the metric under consideration and indices are raised us- 
ing that metric.) For an unperturbed Schwarzschild black hole 
of mass M, PN harmonic coordinates are obtained by transform- 
ing the Schwarzschild radial coordinate V to -Rpn = V — M and 
thus retain the coordinate singularity at the horizon present in 
Schwarzschild coordinates. 
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essary ingredients (except for a few that can be obtained 
from Will and Wiseman [zH). Even more conveniently, 
their expression is in the same (harmonic) coordinate sys- 
tem as BFP's PN metric, so we do not have to determine 
a coordinate transformation for the matching between 
the near and far zones. [We have checked explicitly that 
the near and far zone metrics match through 0(v 5 ) in 
all components.] One also needs to know the binary's 
past history in order to obtain the far zone metric accu- 
rately, due to retardation: We calculate this in the PN 
approximation using Blanchet's results [62| . 

We could have used Alvi's result for the far zone met- 
ric [4l| if we had only been interested in obtaining ini- 
tial data valid through 0(v A ): Alvi's metric, computed 
following Will and Wiseman, is in the same coordinate 
system as Pati and Will's (and thus BFP's near zone 
metric), and is computed through 0(v 5 ) in all its com- 
ponents. However, we needed higher-order results than 
he obtained, or Will and Wiseman's results could give, 
to construct our extended data. This is also why BFP's 
metric was preferable to the ADMTT metric given by 
Jaranowski and Schafer [72j : One needs to calculate sev- 
eral integrals (one of which Kelly et al. [43| had to resort 
to numerical techniques to evaluate) to obtain even the 
0{v A ) pieces of the far zone metric in this approach, while 
it is possible to obtain the 0(v 5 ) pieces of the far zone 
metric merely by taking derivatives using Pati and Will's 
results. [One is also able to obtain initial data through 
0(v 5 ) in the near zone using BFP's results, since they 
give the 0(v 6 ) pieces of the spatiotemporal components 
of the metric. This is not possible with Jaranowski and 
Schafer 's results.] 



2. Comparison with Taylor and Poisson's determination of 
the tidal fields 

Our method for determining the tidal fields can be 
compared and contrasted with that employed in the re- 
cent calculation by Taylor and Poisson [73[ • Most impor- 
tantly, Taylor and Poisson's aims are slightly different 
and more general: They have carried out the matching 
of a single black hole to an arbitrary 1PN metric (ex- 
pressed in terms of potentials), and used the standard 
PN prescription (for obtaining equations of motion) for 
which orders to keep in the various components, as op- 
posed to our "initial data prescription." They thus keep 
terms through 0(v 4 ) in the purely temporal component 
of the metric, 0(v 3 ) in the spatiotemporal components, 
and 0(v 2 ) in the purely spatial components, while we 
keep terms through 0(v 5 ) in all components. [While the 
0(v 5 ) terms in the purely temporal and purely spatial 
components are needed to obtain the 0(v 5 ) piece of the 
coordinate transformation, which itself is necessary for 
obtaining the initial data to a formal accuracy of 0(v 4 ), 
those components themselves do not increase the formal 
accuracy of the resulting initial data.] We also expand 
the near zone metric in multipoles before matching, while 



Taylor and Poisson do not. 

But the general setups have some superficial similari- 
ties: Both approaches use harmonic coordinates for the 
PN metric and THZ coordinates [i^, [7(| for the black 
hole perturbation. However, Taylor and Poisson trans- 
form Poisson's perturbed black hole metric [68[ to the 
THZ gauge instead of using Detweiler's result [5^] di- 
rectly. Additionally, they use PN harmonic coordinates 
for the black hole background, since they have no need 
for their coordinates to be horizon penetrating. 

The actual determination of the tidal fields is car- 
ried out in a very different manner in the two calcula- 
tions: Taylor and Poisson first introduce the most gen- 
eral coordinate transformation that preserves the post- 
Newtonian form of the metric. They then specialize it 
so that it transforms harmonic coordinates to harmonic 
coordinates, and apply it to the PN potentials. After 
this, they decompose the transformed potentials into ir- 
reducible pieces, and can then finally use the matching 
conditions to determine expressions for the tidal fields 
in terms of the potentials. They also obtain the black 
hole's equations of motion and the previously undeter- 
mined pieces of the coordinate transformation from the 
matching conditions. 

We, however, have adopted a more brute force ap- 
proach that assumes nothing about the coordinate trans- 
formation a priori, except its zeroth-order value. (In 
fact, even the zeroth-order value can be shown to be 
constrained by the matching.) However, it does assume 
that the PN point particle trajectories are valid for black 
holes, as is required by the strong equivalence principle: 
Taylor and Poisson demonstrate that this is indeed the 
case (to the order they have matched). Our method also 
requires no decomposition, though it makes heavy use 
of the computer algebra system Maple and the asso- 
ciated tensor manipulation package GRTensopJI [74| . 
We use the gauge invariance of the linearized Riemann 
tensor and linear independence to separate out the por- 
tions of the equations that determine the tidal fields from 
those that determine the coordinate transformation. See 
Sec. IV Al for a detailed presentation of our algorithm. 

With this method, we obtained expressions for the 
tidal fields that can be compared with those of Taylor 
and Poisson: The pieces that we both computed agree 
exactly. See Appcndix lB ll for a comparison, including ex- 
plicit expressions for the tidal fields we obtained. While 
our expressions do not have the full generality of Tay- 
lor and Poisson's, they do include the 1PN corrections 
to the magnetic quadrupole and electric octupole fields 
(for a circular orbit), neither of which Taylor and Poisson 
computed. 



D. Structure of the paper 

We begin by giving an overview of our matching proce- 
dure in Sec. iHland then present expressions for the inner 
and near zone metrics in Sees. IIIII and IIVI In Sec. IIVI wc 
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also consider two relevant aspects of the PN metric, viz., 
conformal flatness breaking and gravitational radiation 
effects. Next we discuss the specifics of our matching 
procedure and read off the matching parameters and co- 
ordinate transformation order-by-order in Sec. |Vj We 
compute the far zone metric in Sec. WI\ where we also 
discuss the PN results we use to obtain the effects of ra- 
diation reaction on the binary's evolution. In Sec. lVIll wc 
give an overview of the construction of an extension of 
this data set that is valid through 0(v 5 ) in the near and 
far zones, in addition to including various other higher- 
order terms. Then we stitch the metrics together numer- 
ically in Sec. IVIIH first resumming the near zone metric 
to improve its strong field behavior (and thus the match- 
ing), then constructing transition functions to stitch all 
the metrics together smoothly, and finally considering 
the constraint violations of the resulting merged metric. 
Lastly, we conclude and summarize in Sec. ILXl 

We present various ancillary results and technical de- 
tails in the appendices: Appendix [2] compares Cook- 
Scheel and PN harmonic coordinates. We provide ex- 
plicit expressions for the tidal fields and some related 
discussion in Appendix [B] along with the calculation of 
the fourth order pieces of the octupole tidal fields and 
the polynomial part of the associated coordinate trans- 
formation. In Appendix [Cl we give the details of our cal- 
culation of the higher-order extension to the data, and 
in Appendix [D] the precise details of how the metrics are 
implemented numerically. 



E. Notation and Conventions 

Units: We use geometrized units with G = c = 1 
throughout. (G is Newton's constant and c is the speed 
of light.) 

Binary parameters: The binary's orbital velocity is v, 
its orbital angular velocity is to, and its coordinate sepa- 
ration is b. The masses of the holes are mi and 7712; their 
total mass is m := mi + m.2- 

Orbit terminology: When we are only considering 
terms through 0{v A ) and the PN equations of motion 
are thus conservative, we shall refer to the binary's or- 
bit as circular. When we are considering nonconservative 
terms [at 0(v 5 ) and higher] but find that their effects can 
be ignored in the current portion of our calculation, we 
shall refer to the binary's orbit as (quasi)circular. When 
we construct the extension to our data in Sec. I VIII and 
Appendix [C] or are considering the binary's past history 
as encoded in the far zone metric in Sec. I VI Al and can no 
longer ignore radiation reaction, we refer to quasicircular 
orbits. 

Indices: We use the standard convention that Greek 
letters denote spacetime indices, while lowercase Roman 
letters (here k, I, p, s, u, and v) denote spatial indices. 
In Sec. IIII1 uppercase Roman letters denote indices on 
the 2-sphcre, but in Sees. HH El and IVIII Bl they label 
the holes, as well as, by extension, the zones into which 



we divide the data's timcslicc. [In the latter context, 
the notation "+ (I <-> 2)" denotes that the preceding ex- 
pression is to be added to itself with the labels I and 
2 switched.] In Sec. I VII and Appendix [Cl Q denotes 
a multi- index. Except where otherwise noted (e.g., in 
Sees. IIV Al IVI[ and IVIII]) . spacetime indices are raised 
and lowered using the Minkowski metric rj a p, so spatial 
indices are raised and lowered using the Kroneckcr delta 
Ski '■= diag(l,l,f) (the symbol ":=" indicates a defini- 
tion). We may even freely raise and lower spatial in- 
dices within expressions for notational convenience, par- 
ticularly in Sec. IIVI The summation convention is al- 
ways in force, and may even be applied to spatial in- 
dices that are at the same level, particularly in Sec. I VII 
Parentheses, square brackets, and angle brackets on in- 
dices denote symmetrization, antisymmetrization, and 
the symmetric trace-free projection, respectively. For in- 
stance, A {af3) := (A a p+Ap a )/2, A [a/3] := {A a p-Ap a )/2, 
and A <u> := {A M + A M )/2 - AP p 5 k i/3. We use verti- 
cal bars to exclude indices from these operations — e.g., 
A( a \p\j) := (Aa^ + A 1 p a )/2. 

Arrays: In addition to the ordinary 3-dimcnsional Kro- 
necker delta defined above, we also define a "lowered 
4-dimensional Kronecker delta," A a p := diag(l, I, I , I). 
Our conventions for the three- and four-dimensional Levi- 
Civita symbols e k i p and e a p 1 s are that €123 = £0123 = 1. 

Metrics: As is usual, r\ a p denotes the Minkowski met- 
ric; our signature is (— , +, +, +)■ In Sees. IIVI and IVl (and 
Appendix [B} g a p denotes the near zone (PN) metric and 
h a /3 denotes the inner zone (perturbed black hole) metric. 
In Sec. lVIIlAl and AppcndixlAl g a p (sometimes with dec- 
orations) is also used for the unperturbed Schwarzschild 
metric. In Sec. IVII (and the associated Appendix [C]) . g a p 
denotes the far zone metric, and h a p its associated metric 
perturbation. 

Coordinates: Inner zone coordinates are X a — 
(T,X k ) : with V := \/X k X k (Schwarzschild); X a = 
(T,A fe ) = (T,X,Y,Z), with R := V ' X k X k (Cook- 
Scheel). PN harmonic coordinates (used in the near 
and far zones) are x a = {t,x k ) = (t,x,y,z), with 
r := \Jx k x k ] r a denotes just the spatial coordinates [i.e., 
r a = (0, x k )]. Unit vectors are denoted by "hats." For 
instance, t a ,x a ,i/ a , and z a are the Cartesian PN coordi- 
nate basis vectors corresponding to indices 0, 1, 2, and 3, 
respectively. Because we are interested in expanding in 
the distance from hole 1, we also define "tilded" coordi- 
nates with their origin at hole l's position at t = 0, viz., 
x a := x a — (m2/m)bx a , x := x — (m,2/m)b, r := y x k x k , 
and f a := r a — (rri2/m)bx a . Spatial vectors (or the spa- 
tial parts of spacetime vectors) will be denoted cither 
with an arrow or a spatial index. 

Norms: The Euclidean norm for spatial vectors is de- 
noted by || • || (so, e.g., we could write the definition of R 
above as R :— \\X\\). 

Derivatives: All partial derivatives are taken with re- 
spect to harmonic PN coordinates, so 9„ := d/dx a . 
Overdots denote differentiation with respect to t (i.e., 
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PN harmonic time). While in Sec. IIII1 ovcrdots properly 
denote differentiation with respect to T (i.e., the Cook- 
Scheel time coordinate), as we shall see, this is equivalent 
to differentiation with respect to t to the order we are 
considering. We use the shorthand d a p := d a dp. 

Order counting: We have to deal with two expansions 
here, since we are performing asymptotic matching. For 
our purposes, the PN expansion can be treated as an 
expansion in y/ m^/b. We use m,i (instead of m,\ or m) 
for convenience, as we concentrate on matching around 
hole 1. (We can do this without loss of generality, by 
the symmetry of exchanging labels 1 and 2.) The black 
hole perturbation expansion is in multipoles, so around 
hole 1 it can be treated as an expansion in r/b. Here we 
treat t/b as the same order as f/b, as required by the slow 
time variation assumption made in the derivation of the 
perturbed black hole metric. This means that we only 
match in a neighborhood of t — 0. 

Order notation: It will be important for us to have a 
compact notation for, e.g., the jth term of the expansion 
of some arbitrary quantity. For this, we use the "order 
projection" notation, which has three different flavors: 
Since the slow-motion PN expansion can be considered 
our primary expansion (and is the expansion to which 
we refer when we simply call something jth order), the 
most commonly used flavor will be (•) .. This denotes the 

coefficient of {m^/by^ 2 in the (asymptotic) power series 
expansion of its argument. For the reasons discussed in 
Sec. IV Al this includes multipoles through octupolar or- 
der for j 6 {2, 3} but only through quadrupole order for 
j € {4, 5}. (In principle, we include all the multipoles for 
j € {0, 1}, since there are only monopole contributions.) 
For the times when we need to single out a particular mul- 
tipole for consideration, we shall use (•) ■ n , which denotes 
the 2™th multipole-order piece of (•) .. Finally, when we 
want to include all of the multipolar orders through 2" 
we shall use (-)j,< ra - 

Polynomial and nonpolynomial parts: All the expres- 
sions we consider when performing the matching of the 
inner and near zones can be written as the sum of a 
polynomial in x a and nonpolynomial terms of a par- 
ticular form, viz., polynomials in x a multiplied by f™, 
n € Z\ {0, 2,4,.. .}. It will often be convenient to divide 
expressions up into their polynomial and nonpolynomial 
parts, since these can be treated separately, by linear 
independence. We shall thus denote the polynomial and 
nonpolynomial parts of an expression by the superscripts 
P and NP, respectively. 

References: (3?| and [H[ will be referred to as Papers I 
and II, respectively. 

II. ASYMPTOTIC MATCHING 

The technique of asymptotic matching is a standard 
one in the analysis of multiscale and singular perturba- 
tion problems, allowing one to relate and combine ap- 
proximate solutions that are valid on different scales [75| . 



It has been used in general relativity to obtain, e.g., PN 
equations of motion — [73| contains the most recent of 
these calculations — and the radiation zone metric of a 
binary system (76[. (See Paper I and [73[ for further dis- 
cussion and references.) Here we specialize our discussion 
to the case of a binary black hole. 

A timeslicc of a binary black hole spacetime divides 
naturally into four primary zones and three secondary 
buffer zones: See Fig. [5] for an illustration. In practi- 
cal work, the boundaries of all of these zones are nec- 
essarily only given approximately, since we do not cur- 
rently possess sharp estimates for the error bars of the 
approximations used to describe this spacetime. There 
are two inner zones around the black holes, given by 
ta <C b, where is the distance from (the point parti- 
cle associated with) hole A: Here the spacetime is well 
described by a perturbed black hole. Surrounding (and 
partially overlapping) these is the near zone, where the 
standard (harmonic coordinate) PN metric is valid — i.e., 
not too close to the holes, yet not so far away that retar- 
dation cannot be treated pcrturbativcly. This is given by 
rA 3> rriA and r < A, where tua is the mass of hole A, 
r is the distance from the system's center-of-mass, and 
A = b/2v is the reduced characteristic wavelength of the 
binary's gravitational radiation. Finally, the remainder 
of the timeslice (along with a little bit of the outer por- 
tion of the near zone) comprise the far zone, given by 
r > A, where retardation can no longer be treated per- 
turbativcly, and spacetime is described by a separate PN 
metric that accounts for this. (N.B.: Different relations 
were given for the outer edge of the near zone and inner 
edge of the far zone in Papers I and II, but the ones used 
here are more accurate.) See Table IIIII (in Sec. IVIII Bp 
for numerical values for the boundaries of the zones for 
a particular equal-mass binary. 

The three buffer zones are the portions of the times- 
lice where the preceding four regions overlap. (Due to 
the "fuzziness" inherent in < and >, the near and far 
zone can have a substantial overlap, despite formal ap- 
pearances.) We restrict our attention to cases in which 
the specified zones overlap in the manner we shall de- 
scribe, but in no other fashion (e.g., we do not want the 
two inner zones to overlap). There are two buffer zones 
where the two inner zones overlap the near zone (thus 
these are given by niA <C &), and another buffer 

zone where the near zone overlaps the far zone, which is 
very roughly an annulus whose radius and thickness are 
both of order A. (The order of the thickness of the an- 
nulus is our choice — see Sec. IVIII Bl All that is required 
formally is that it increase as v decreases, correspond- 
ing to a larger realm of validity of the near zone metric's 
pcrtubative treatment of retardation. Additionally, as 
is indicated in the figure, this buffer zone would not be 
spherical in a more nuanced description.) These buffer 
zones are where we perform the formal matching that 
determines the coordinate transformation and relations 
between parameters as well as where we stitch the met- 
rics together numerically. 
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FIG. 2: A diagram of the zones into which we divide the 
binary black hole timeslice. The two black holes (BH 1 and 
BH 2, denoted by filled-in circles) lie on the x-axis, surrounded 
by their respective inner zones (Ci and C2) and inner-to-near 
buffer zones (O13 and 023)- (In actuality, the black holes 
should be tidally distorted, along with their associated inner 
and buffer zones. We neglect this distortion in the above 
diagram for simplicity.) The near zone, C3, covers the orbit 
of the binary and is surrounded by the far zone, C4, and the 
near-to-far buffer zone, O34. 



This matching and stitching relies on the observation 
that if the metrics that are valid in the various zones are 
all to be different approximations to the same (unknown) 
global exact metric, then, considered as abstract tensors, 
they should agree with each other in their realms of mu- 
tual validity. More specifically, assume that there exists 
a buffer zone in which two approximate metrics g' 1 ' and 
g( 2 ) are both valid. (We write these metrics without 
indices to emphasize that they are currently being con- 
sidered as abstract tensors.) Take their associated small 
parameters to be t\ and 62, respectively. (The buffer zone 
is then defined to be the region in which these parame- 
ters are indeed small.) Then make a bivariate expansion 
of the metrics in both small parameters. That is, take 
gW, which is already an expansion in ex, and expand it 
in £2 as well; similarly expand g*- 2 ) in e\. The coefficients 
of both bivariate expansions, considered as abstract ten- 
sors, should then be equal if the metrics describe the same 
spacetime. The resulting equations relate the parameters 
of the two metrics. 

While the statement of this result in terms of abstract 
metrics is simple, in practice one works with the metrics' 
coordinate components. In general, the coordinate sys- 
tems in which the metrics' components are known will 
not agree, so one first chooses the coordinate system in 
which one of the metrics is expressed to be the primary 
coordinate system. (For us, this will be the near zone's 
PN harmonic coordinate system.) One then applies an 
arbitrary coordinate transformation to the other metric 
in order to put it in the same coordinate system as the 
first (to the order one matches). Equating coefficients of 
the bivariate expansions of the components of the metrics 
(including the coordinate transformation) is then equiv- 



alent to equating the coefficients of the expansions of the 
abstract tensors. In this case, the resulting equations 
will determine not only the relations between the met- 
rics' parameters, but also the arbitrary coordinate trans- 
formation. In our calculation, each order's contribution 
to the coordinate transformation contains some arbitrary 
pieces — i.e., pieces that can be set to any value without 
affecting that order's matching. We have found that most 
of these are determined by the matching at two orders 
higher. 



III. INNER ZONE METRIC 

Detweiler's perturbed black hole metric [59| is derived 
under the slow-motion and weak-curvature assumptions 
of Thorne and Hartle [6^|. Here one characterizes the 
strength of the perturbation by the length scale 7Z, de- 
fined to be the smallest of the radius of curvature, in- 
homogeneity scale, and scale of time variation of the ex- 
ternal universe. For convenience, these are generally all 
taken to be formally of the same order. (In the nomencla- 
ture of Thorne and Hartle, the "external universe" refers 
to the perturbing spacetime in which the black hole is 
placed, here consisting of the tidal fields generated by the 
hole's binary companion.) For a post-Newtonian binary 
we have 1Z ~ \/h z Im-i [see, e.g., the discussion following 
Eq. (1.4b) in [69||]. Here b is the coordinate separation 
of the holes and 7712 is the mass of the hole's companion 
(the one responsible for the perturbing tidal fields). (We 
are thus specializing to the inner zone surrounding hole 1 
in giving this scaling, as well as the scalings of the tidal 
fields below. One can obtain the scalings for the inner 
zone surrounding hole 2 using the substitution 1 <-> 2.) 

With these assumptions, the time derivative of a 
quadrupolc field is of octupolar order. Since one is only 
able to include quadrupole and octupole perturbations 
in linear black hole perturbation theory, the time depen- 
dence of the tidal fields in the metric is restricted to the 
quadrupole fields, and even there it is only linear. Of 
course, it is possible to obtain the complete time depen- 
dence of the tidal fields (for times much less than the 
radiation reaction timescale), as Taylor and Poisson [73[ 
did, since we can compute the tidal fields at any point in 
the orbit. (See Appendix [Bl for further discussion.) How- 
ever, knowing this dependence does not let us compute 
the metric to a higher formal accuracy. 

In Detweiler's metric, the tidal perturbations are en- 
coded by symmetric trace-free electric and magnetic 
tidal fields £m (electric quadrupole), £ki p (electric oc- 
tupole), Bki (magnetic quadrupole), and Bki P (magnetic 
octupole). These come from the Thornc-Hartle-Zhang 
(THZ) harmonic specialization of locally inertial coor- 
dinates (from [gg, [z3|), which Detweiler uses to express 
the perturbation of Minkowski space that his metric ap- 
proaches in the buffer zone. 

As their names suggest, the quadrupole and octupole 
tidal fields formally scale as TZ~ 2 and 1Z~ 3 , respectively. 
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However, this is only a formal scaling, since it requires 
that one treat the radius of curvature, inhomogeneity 
scale, and scale of time variation of the external uni- 
verse as being of the same order [see, e.g., Eq. (1.4b) 
in [69(| ] . For a post-Newtonian binary, these are not 
of the same order: The magnetic moments are a fac- 
tor of v smaller than the corresponding electric mo- 
ments; the octupole moments are a factor of v larger 
than the formal scaling would suggest — see Eqs. (2.1) 
and (2.2) in [73[ and the surrounding discussion. Specif- 
ically, the magnetic quadrupole tidal field on hole 1 is 
actually ~ (m 2 /b) 1 / 2 (m 2 /b 3 ) [cf. the expression given in 
Eqs. (5.45b) and (5.56b) of [23], and the electric octupole 
tidal field on hole 1 is actually ~ m 2 /b 4 [cf. Eqs. (5.50) 
and (5.56a) in [z3|]- In order to simplify later expressions 
and make their expansion in r/b explicit (f denotes the 



distance from hole 1), we define "barred" versions of the 
tidal fields around hole 1 with the scaling taken out, viz., 

Sh := (b 3 /m 2 )£ kl , B kl := (b/m 2 )^ 2 (b 3 /m 2 )B k i, 

Skip ■= (b 4 /m 2 )£ki p , Bkip := {b/m 2 ) 1/2 (b 4 /m 2 )B k i p . 

(3.1) 

Detweiler (59j presents the metric perturbation in 
Schwarzschild coordinates, giving the portion involving 
the quadrupole fields, including their time derivatives 
(which are actually of octupolar order), in his Eqs. (G.6)- 
(G.ll), and the portion with the octupole fields in his 
Eq. (58). Putting this all together, along with the back- 
ground unperturbed Schwarzschild metric, we obtain a 
line element of 



h i ^dX a dX f3 = -H^dT 2 + H ( k S) dTdX k - - 
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(3-2) 



where the "(S)" acts as a reminder that this is in Schwarzschild coordinates X a = (T, X k ), with a radial coordinate of 
V := ^X k X k . Additionally M is the mass of the hole, (tab denotes the metric on the 2-spherc, and we have defined 
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(3.3) 



We have also defined C k i p :— e k i s B s p and C k i ps := 
e k i u B u ps , for convenience, as the magnetic tidal fields 
only appear in the perturbation in these dual forms. (The 
"barred" versions of these quantities are defined analo- 
gously to B k i and B k i p .) Here e k i p is the 3-dimensional 
Levi-Civita symbol, with 6123 = 1. While we have not 



indicated this explicitly, £ k \ and C k i p both depend (lin- 
early) on the null ingoing Eddington-Finkelstein coordi- 
nate V := T + V + 2M\og(P/2M - 1). (The octupole 
fields are treated as constants, since their time deriva- 
tives are of hexadecapole order.) Thus, primes denote 
derivatives with respect to V. We include the B' u contri- 
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butions here, even though they do not increase the for- 
mal accuracy of our initial data: We fix the lowest-order 
piece of B' kl when we read off the 1PN correction to the 
electric octupole (in Appendix IB 2|) as an application of 
our matching procedure. We also are able to obtain the 
full time dependence of the tidal fields (up to radiation 
reaction effects) a posteriori — see the discussion in Ap- 
pendix IB 21 (Including such formally higher-order pieces 
can improve the accuracy of the data in practice.) 

There is one gauge subtlety here that Detweiler [5§| 
does not address in his paper: If one compares the ex- 
pression for the THZ tidal perturbation in his Eq. (53) 
to that of the metric perturbation in his Eqs. (56)-(57), 
it appears that the portions of the metric perturbation 
involving time derivatives of the tidal fields do not ap- 
proach the analogous portions of the THZ tidal pertur- 
bation. The resolution of this apparent discrepancy is 
that Detweiler made a tacit gauge transformation (away 
from the pure THZ gauge, though only affecting the time 
derivatives of the tidal fields, to the order considered) to 
obtain the given compact expressions for the metric per- 



turbation [781 ]. 

We now transform this metric to the quasi-Cartesian 
form of Cook-Scheel harmonic coordinates (from [671] ). 
which we denote by X a = (T,X k ), with R := \fX k X k . 
The transformation is given by 



T =T — 2M log 



R - M 



R+M 



X 



M 
1+ R 



(3.4) 



X k 



so V = R + M. This comes from Cook and Scheel's 
Eqs. (20), (41), and (43) upon noting that Boyer- 
Lindquist coordinates reduce to Schwarzschild coordi- 
nates for a Schwarzschild hole. In order to simplify our 
transformation, we use the relation V 2 aAsdX A dX B = 
(1 + M/R) 2 (dX p dXP - dR 2 ), obtained from dX k dX k = 
dV 2 + P 2 aAsdX A dX B , to eliminate cab, thus obtaining 
a full transformed line clement of 
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(3.6) 



3 R + M, 

We have used the fact that the time dependence of £ k i and C k i p is only linear (to the multipolar order we are 
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considering) to write £' M = £m and £ki(V) = £ u + [T + R + M + 2M \og(V /2M)]£ki, with analogous expressions 
involving Cki p - [This expansion removes the logarithms present in the expression for the metric in Eq. (|3.2p .] Here 
an overdot denotes a derivative with respect to T, and the tidal fields on the right are constant, so £ki and £ki are 
treated formally as independent tidal fields, despite the notation. That is, £ki just denotes the constant part of 
£m{T) = £ki + T£]~i. Finally, despite appearances, this expression for the metric is in fact in a quasi-Cartesian form: 
dR should just be considered a shorthand for X^dX k / R. 



IV. NEAR ZONE METRIC 

We take the harmonic coordinate metric to the order needed from Eqs. (7.2) in Blanchet, Faye, and Ponsot [55 
and specialize it to a circular orbit, obtaining 
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Here xa(0i ^ 6 {1, 2} denotes the position of (the point 
particle associated with) hole A. 3 Thus ?a :=- x — XA(t) 
gives the displacement from hole A, with ta '■= 
giving the distance from hole A and fiA '■= TaI^a the 
associated unit vector. Similarly, va '■= xa denotes the 
velocity of hole A. The displacement vector from hole B 
to hole A is given by tab '■= Ta — ?b, with an associ- 
ated unit vector of tiab '■= Tab/\\^Ab\\- Thus, specializ- 
ing to a (quasi)circular orbit, ||ri 2 || = b, where b is the 
separation of the holes and is therefore constant up to 
orbital shrinkage (which we can neglect when performing 
the matching, given the order to which we are calculat- 
ing, though we include it when implementing the metrics 
numerically — see Sec. IVI A|) . Additionally, we shall usu- 
ally use the shorthand b := f\ 2 . (We do not use it in the 
above expression for the metric since r*i 2 changes sign un- 
der 1 «-> 2, while b does not.) Similarly, vab '■= va — vb; 
we also have S := n + r 2 + b. Finally, the notation 



3 We shall henceforth refrain from belaboring the distinction be- 
tween the PN point particles and true black holes, except where 
we feel that it is important to emphasize this point. 



"+ (1 <-> 2)" denotes that the preceding expression is to 
be added to itself with the labels 1 and 2 switched. 

For a circular orbit, we have b = —uj 2 b, by definition, 
where 

,2 1 



u> 



m 
1? 



(4.2) 



is the (harmonic coordinate) angular velocity [obtained 
from Eq. (8.6) in [58||], and r\ := toito 2 /to 2 is the sym- 
metric mass ratio. Assuming a (quasi) circular orbit [so 
that the separation vector of the point particles is orthog- 
onal to their velocities up to 0(v 5 ) radiation reaction ef- 
fects] also allows us to make the simplifications ni 2 • vi = 
h\ 2 -v 2 = h\ 2 -v\ 2 = 0(v 5 ) in obtaining the expression for 
the metric given in Eqs. (|4.ip . We take the explicit ex- 
pression for the orbit to be b = b(x cos ujt+y sin cut). From 
the definition of the center-of-mass coordinates used in 
the PN metric, the positions of (the point particles as- 
sociated with) the holes are given, to the order we need 
them, by 
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[For a circular orbit, the first PN corrections to these are 
0(t! 4 ) and are thus not needed here.] The expressions 
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for the holes' velocities are obtained by taking a time 
derivative. 

The explicit expansions of the near zone metric at var- 
ious orders given in the next section were obtained using 
the expansions 
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Here we include the octupole pieces in the expression 
for 1 /V2 because they become of quadrupole order when, 
e.g., multiplied by m\jr\ in the — 2miTO2/rir2 piece of 
goo ■ We do not need to expand 1/S n to octupolar order 
as it is not multiplied by mi/n at the PN order we are 
working. 



A. The PN metric and conformal flatness 

At 0(v 4 ), the spatial PN metric is no longer confor- 
mally flat: That is, there does not exist a coordinate sys- 
tem in which the 3-metric can be written as g k i = ^S^i- 
This can be guessed from a cursory inspection of the 
PN metric, as nondiagonal spatial components first ap- 
pear at 0(v 4 ), but was established more firmly in [73 |. 
and then further studied in [soj . In both cases, this re- 
sult was an offshoot of a comparison of the predictions 
of a post-Newtonian analysis of the Isenberg-Mathews- 
Wilson (IMW) approximation [U H3] with those of the 
full PN approximation: They were found to differ start- 
ing at 0(i> 4 ). Such a comparison gives the desired result 
because the IMW approximation assumes spatial confor- 
mal flatness (along with maximal slicing — i.e., a vanish- 
ing trace of the extrinsic curvature) in an attempt to 
remove the dynamical degrees of freedom from the grav- 
itational field. 

However, the comparison with the IMW approxima- 
tion is not exact. The post-Newtonian analyses use the 
ADMTT gauge, for which we have 8ki^ kl = 0, where tt m 
is the canonical momentum conjugate to the 3-metric, 



c/fcz, while maximal slicing sets g lp Ki p = 0, where Ki p is 
the extrinsic curvature of the slice. Of course, these two 
slicings agree as long as gu is a multiple of the flat space 
metric, since then 6kiir kl — gu^ kl = g lp K lp = 0. 
[See, e.g., Eq. (E.2.31) in Wald @ for the relation be- 
tween ir lp and K lp ] But in general, the two slicings will 
differ at the higher, non-conformally flat orders we are 
interested in. 

It is also possible to demonstrate this lack of spatial 
conformal flatness directly, and we shall do so here for 
the harmonic slicing we are using. In four or more di- 
mensions, the Weyl tensor settles questions of conformal 
flatness: It vanishes if and only if a space is conformally 
flat (84[. However, in three dimensions the Weyl tensor 
vanishes identically, and its analogue for settling ques- 
tions of conformal flatness is the Bach or Cotton- York 
tensor, Cki- This is defined (with indices raised by the 
3-metric) by 

C H := 2e k ps V s R lp - ^e a p V p R. (4.6) 

Here Vfc, Rm, and R are, respectively, the (three- 
dimensional) covariant derivative, Ricci tensor, and Ricci 
scalar associated with g^i- In fact, the nonvanishing of 
the Cotton- York tensor is a necessary and sufficient con- 
dition to render its associated 3-metric non-conformally 

flat ( for a P r00I > sec Cha P- VI i § 5 in HI)- 

We have computed the formal 0(c~ 4 ) portions of this 
tensor [i.e., those that correspond to the 0(c -4 ) 0{v A ) 
pieces of the spatial metric] symbolically using Maple 
and GRTensorII and verified that certain components 
arc nonvanishing at various points of the timeslice. As 
an illustration, we have plotted the lowest-order piece 
of the norm of the Cotton- York tensor, ^ 8 k v8 ls CkiC ps 
along the axis passing through the holes in Fig. [3] for the 
standard equal-mass test system of mj = m 2 = m/2, 
b = 10m. As expected, the values are largest in the 
region around the holes, showing that this is the region 
in which the largest perturbation would be required to 
make the 2PN metric conformally flat. 

There is, however, a more intuitive way of understand- 
ing the breaking of spatial conformal flatness of the PN 
metric (here returning to the ADMTT slicing). This ap- 
proach relates explicitly to the failure of the PN metric to 
be manifestly conformally flat at 0(v 4 ) and comes from 
the work of Nissanke [42j . 4 From the study of PN theory 
in the canonical ADM formalism [87| , one knows that the 
PN metric can be rewritten as 



9ki 
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4 This argument was suggested to one of us by Luc Blanchet l86l . 
Additional intuition can be obtained from Valiente Kroon's 
result [4(ill that stationary spacetimes only admit conformally 
flat slices if certain "obstructions," constructed out of Geroch- 
Hanscn quadrupolcs, vanish. 
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FIG. 3: The norm of the Cotton- York tensor for the test 
system mi = — m/2, b — 10m along the a;-axis (i.e., the 
axis passing through both holes). (We only show it around 
hole 1 because it is symmetric about x = 0.) 



where 'F is some conformal factor, while hki is not propor- 
tional to the flat metric. In fact, hki must be symmetric 
and trace-free and contains a piece that is proportional 
to 



dx< 



<k dx l > 



log 5, 



(4.8) 



where S was defined in the previous subsection and the 
angle brackets stand for the symmetric trace-free projec- 
tion. [See, e.g., Eq. (3.1) in [42|]; we have left off pieces of 
the form v"^ k v l {* .] This is clearly not a manifestly con- 
formally flat contribution to the metric, so we can use 
it to connect with the heuristic that "the 2PN metric is 
not spatially conformally flat because it contains various 
pieces that are not manifestly conformally flat." 

To do so, we shall demonstrate how the presence in 
hki of the term given in Eq. (|4.8[) prevents the obvi- 
ous sort of coordinate transformation from rendering 
conformally flat. For the purposes of illustration, we 
take our coordinate transformation to be of the form 
x k — > x' k = x k + where £ fc = 0(v 4 ), and ask that 
it remove the h^i piece from the 3-metric. (Of course, in 
general the coordinate transformation would just need to 
turn hki into a scalar multiple of Ski , and would not have 
to be of the form above, but our assumptions suffice for 
a heuristic argument.) Thus, £ fe should satisfy 
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but this equation has no solution when the right-hand 
side is given by Eq. (|4.8[) . To show that this is the case, it 
is sufficient [by the gauge invariance of the linearized Ric- 
mann tensor — see the discussion surrounding Eq. (|5.5|) ] 
to show that the (three-dimensional) flat space linearized 
Riemann tensor associated with hki does not vanish. We 
have computed this tensor [with hki given by Eq. (|4.8[) ] 
and verified that several components are indeed gencri- 
cally nonzero. 



B. Gravitational radiation in the PN metric 

The order at which the effects of gravitational radi- 
ation appear in the PN metric is different for different 
effects, which can easily lead to confusion. We thus offer 
a brief discussion of these orders here. To avoid even fur- 
ther confusion, since we shall discuss some pieces that are 
not dimensionless, we shall describe all orders in terms of 
the formal slow motion expansion in 1/c, as do Blanchct, 
Faye, and Ponsot (BFP) (58j[. [Since the metric is di- 
mensionless, an 0(c~ n ) contribution to it can be unam- 
biguously identified as being 0(v n ).] It is well known 
that the effects of gravitational radiation reaction first 
enter the equations of motion at Q(c~ 5 ) for a circular 
orbit. [See, e.g., Eq. (189) in Blanchet [62||.] This is also 
the leading order of the binary's gravitational wave lu- 
minosity. [See, e.g., Eq. (171) in Blanchet.] However, 
the lowest-order "quadrupolc formula" piece of the grav- 
itational waveform appears in the PN far zone metric 
at one order lower, viz., 0(c~ 4 ). This can be seen in 
Blanchet's Eq. (238); his x variable is 0(c~ 2 ). These 
terms can also be seen in Eqs. (6.10) and (6.11a) of [7l| . 
Here they are presented in a form that allows for a more 
direct comparison with our expression for the far zone 
metric, though this expression does not show the factors 
of c _1 explicitly. The terms in question appear as the 
final two terms in the curly brackets in our expression 
for the far zone spatial metric in Eq. (|6.6c[) . Explicitly, 
they are 4(mi/r)(wf - u> 2 x kt ) + (1 <-> 2). 

Our initial data should thus contain the binary's out- 
going gravitational radiation (whether one uses the ex- 
tension that adds on various higher-order terms or not), 
since we have included all the 0(v 4 ) pieces in the spatial 
metric, along with the matching 0(v 5 ) pieces in the ex- 
trinsic curvature. This explains why Kelly et al. obtain 
the outg oing wave train when they evolve [571 ] , since their 
data set [431 ] contains only terms in the far zone through 
0{v A ) [in the spatial metric, with the matching 0(v 5 ) 
pieces in the extrinsic curvature], even though they do 
not have the 0(v 5 ) terms in the spatial metric. 

The factor of 1 /c one obtains in going from the terms 
in the near zone that give the lowest-order piece of 
the waveform when expanded in the far zone (i.e., in 
b/r), which are 0(c~ 4 ), to the lowest-order 0(c~ 5 ) ra- 
diation reaction contribution to the equations of motion 
is explained by Blanchet in [88| : Radiation reaction ef- 
fects arise from antisymmetric waves — i.e., those that are 
time-odd — and such waves will come from expressions in- 
volving at least one more time derivative (and thus factor 
of 1/c) than their symmetric (and time-even) counter- 
parts. Explicitly, these near zone terms that match onto 
the lowest-order gravitational wave terms in the far zone 
are A(mi / r\)v k v[ — 4(mirri2 / bS)n k 2 n 12 + (1 * _> 2), since 
r A = r[l+0(b/r)} and S = 2r[l+0(b/r)]. The (multipo- 
lar expansion of the) first time derivative of these terms 
gives the lowest-order radiation reaction contribution to 
the metric given in [88| . These terms themselves are also 
some of the contributions to the near zone metric that 
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cause it to be conformally curved. 

V. THE MATCHING CALCULATION 
A. The setup 

By symmetry, we can concentrate on performing the 
matching around hole 1: All our results in this case can 
be directly translated to the matching around hole 2 
by taking 1 <-> 2, (t,x,y,z) — > (t,—x,—y,z), and 
(T,X,Y,Z) -> (T, -X, — Y, Z), along with x a -> -i a 
and y Q — > — y a in the tidal fields. Here x a and denote 
the unit vectors in the x\- and ^-directions, respectively. 
(Similarly, t a and z a denote unit vectors in the xq- and 
j^-directions — we shall use these later.) In other words, 
we switch the masses to turn mi into rri2, and then ro- 
tate by 7r radians around the z-axis to move the new m-i 
into the same position as the old ra^- However, we also 
need to rotate the inner zone coordinates so that they 
have the same relation to the new near zone coordinates 
that they did when we performed the original match- 
ing around hole 1. Finally, the tidal fields transform as 
Cartesian tensors under this rotation, and this transfor- 
mation is taken care of by the final two substitutions. 

Since our matching calculation will determine the coor- 
dinate transformation, the relations between the metrics' 
mass parameters, and the inner zone mctric's tidal fields, 
we need to posit expansions for all of these. For the coor- 
dinate transformation, we make nearly the same ansatz 
as in Papers I and II, viz., 

X <V) = J2 (^) (*"), + 0(v% (5.1) 

(This is slightly more general than the ansatz used in the 
previous papers because we do not fix the zeroth order 
coordinate transformation from the outset.) We choose 
yrwz/b as our expansion parameter because it makes 
for slightly simpler notation than either v = drnjb or 
■sj mi/b when expanding in the buffer zone around hole 1 , 
as we are doing here. We can make this choice without 
loss of generality: The resulting coordinate transforma- 
tion will be the same regardless of which of these three 
possibilities we choose to use as our expansion parameter, 
though the coefficients of the expansion parameter will 
differ by ratios of the masses. As in the previous papers, 
we implicitly assume that (X a )j is a power series in f/b, 
including negative powers — i.e., a Laurent series — so we 
can write, e.g., 7712/r = (rri2/b)(b/r). However, {X a )j 
should not depend on rri2/b, by definition. 

We also need to worry about the multipole expansion 
of each {X a )-. This would seem to be straightforward, 
since we only want to keep terms through quadrupole 
order overall. However, the structure of the inner zone 
metric creates some complications: In order to obtain 
data that include all the quadrupole [0([f/&] 2 )] pieces 



at fourth and fifth orders, 5 one needs to obtain the oc- 
tupole [0([^/^] 3 )] pieces of the coordinate transforma- 
tion when matching at second and third orders. This is 
due to the appearance of b/f terms in (h a p). for j > 2. 
These enter the fourth and fifth order coordinate trans- 
formation equations, where they multiply the second and 
third order pieces of the coordinate transformation and 
thus produce quadrupolar contributions from octupolar 
pieces of the coordinate transformation. The octupolc 
fields themselves also enter, as they are multiplied by 
b/f in the fourth and fifth order pieces of the inner zone 
metric. However, we shall see that the octupole piece 
of the second order coordinate transformation vanishes, 
and that of the third order one only appears in the time 
coordinate. Thus, the only place where the octupolc 
pieces of the coordinate transformation appear is in the 
(/i( Q | 7 ) 2 (^4|^) 7 ) 3 piece of the fifth order coordinate trans- 
formation equations [given in full in Eq. (|5.30l) ]. 

This increase in the number of multipoles that have to 
be kept as one proceeds to higher and higher orders in 
v is a general feature of the matching of these two met- 
rics. It is thus a source of significant technical difficulty: 
One would need to include the hexadecapole pieces in 
the matching calculation if one wanted to include all of 
the quadrupole pieces at sixth and higher orders. This 
follows since at sixth order the hexadecapole pieces start 
to be multiplied by 6 2 /f 2 , making them of quadrupole 
order. (It would still be possible to obtain all the dipolc 
pieces at fifth and sixth orders if one had an expression 
for the near zone metric through seventh order, but one 
would not be able to include any new higher-order cor- 
rections to the tidal fields this way.) 



We posit the same expansion in \fm2Jb for the mass 
parameter of the inner zone metric, Mi , as we did for the 
coordinate transformation, so 



M 1 =^(^y /2 (M 1 ) i +0(A (5.2) 

J=0 



However, it will turn out that we did not need to allow 
this freedom, as we shall find that the mass parameters 
of the two metrics agree to the highest order to which our 
matching fixes them — i.e., Mi = mi + 0(v 4 ). Similarly, 



We count orders using our primary expansion in v (or, equiva- 
lently, *Jma/V). 
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we asymptotically expand the tidal fields in \J m^Jb, so 

3=0 

3=0 

3=0 

(5.3) 

[The expansions of the duals of the magnetic fields — 
e.g., Chip — arc defined analogously. Also recall that the 
"barred" tidal fields are defined in Eq. (|3.1|) . Addition- 
ally, the tidal fields and their time derivatives — e.g., Ski 
and Ski — are treated as formally independent.] We do 
not include the 0(v 5 ) pieces of the octupolc tidal fields 
here because we did not fix them in the matching: We 
only had to match the octupole fields through 0(v 3 ) to 
obtain initial data with formal uncontrolled remainders 
of 0(v 5 ) and 0([r/b} 3 ) (i.e., octupolar order). We also 
chose to read off the 0{v A ) pieces of the octupole fields 
separately (in Appendix IB 2[) . but did not do so for the 
0(v 5 ) parts. 

In order to read off the matching parameters (and any 
undetermined pieces of the lower-order coordinate trans- 
formations) as efficiently as possible, we note that all of 
our equations for the coordinate transformation at orders 
beyond the zeroth will be of the form 

d{ a Xp) = S aj 3 = S( a py (5.4) 

Here d a := d/dx a (i.e., all partial derivatives are taken 
with respect to PN harmonic coordinates), X a is a func- 
tion of x a , and S a p, the equation's source, is some (sym- 
metric) matrix function of x a (either explicitly, or implic- 
itly though X a ) which is C 2 in the buffer zone. The inte- 
grability condition for this equation is that the flat-space 
linearized Ricmann tensor associated with S a p vanish, 
i.e., that we have 

1 a f3~j& '■— dapS^s + d 1 gS a fi — d a sS 1 (3 — djfjSas — 0. (5.5) 

[N.B.: For convenience, we have defined lap-yS with a 
slightly different index ordering than the linearized Ric- 
mann tensor — given in, e.g., Eq. (5.44) of [77| — and with- 
out the factor of 1/2.] This follows from the gauge invari- 
ance of the linearized Ricmann tensor. (See, e.g., §4.1 in 
Straumann [89| for a proof and discussion of that result. 



In [9(j, Blanchet and Damour use this gauge invariance 
for the same purpose we do.) 

For future use, we note that the homogeneous equa- 
tion, d( a Xp\ — 0, is the flat space Killing equation, and 
its most general solution is given by 

X a = F Pa x p + C a , (5.6) 

where F a p = -F[ Q| gi is some constant (antisymmetric) 4x4 
matrix and C a is some constant 4x1 matrix. See, e.g., 
§13.1 in Weinberg [9l[ for a proof. (Our X a , F a p, and 
C a correspond to Weinberg's £ a , bp a = —b a p, and a a , re- 
spectively.) We shall primarily employ this result tacitly 
at each order beyond the zeroth to ensure that we have 
the most general expression for that order's contribution 
to the coordinate transformation. 

Our general approach to the nontrivial matching that 
occurs at second order and beyond will be as follows: 
We first use the above integrability condition to read off 
the matching parameters, exploiting the linear indepen- 
dence of various terms to simplify the process and justify 
our claims of uniqueness. We start with the nonpolyno- 
mial terms, which determine many of the previously un- 
determined parts of the coordinate transformation from 
two orders lower, as well as the inner zone mass parame- 
ter; the polynomial part then determines the tidal fields. 6 
The nonpolynomial part consists of all the terms that are 
not polynomials in x a := x a — (m2/m)bx a . For this cal- 
culation, these arc all of the form of a polynomial in x a 
multiplied by f n , where n G Z \ {0, 2, 4, 6, • ■ ■ }. After we 
have fixed all the parameters that can be fixed at a given 
order, we can then solve for that order's contribution to 
the coordinate transformation (and the polynomial part 
can be solved for separately from the nonpolynomial part 
that first appears at fourth order). In all of this, Maple 
and GRTensorII proved very helpful: They work ex- 
tremely well for all aspects of the polynomial part, while 
requiring more care when applied to the nonpolynomial 
part. 



B. Zeroth Order [0([m 2 /b] )] 
To lowest (zeroth) order in (n^/fr) 1 / 2 , we have 

(ffa/Oo = (A* 7 ) (Mo (V)„ > ( 5J ) 

where we have defined A a @ := d a X@. Since (g Q ^) = 
(h a /3)o = Va/3> matching just tells us that (A a /3 ) Q must be 
a (general) Lorentz transformation — i.e., not necessarily 



It is intuitively reasonable that the polynomial and nonpolyno- 
mial parts should determine the parameters that they do: If 
one neglects all gauge subtleties and the like, the nonpolynomial 
terms can be thought of as being those associated with hole 1, 
and the polynomial terms with (the tidal fields of) hole 2. We 
discuss this more fully in Sec. IV Gl 
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one continuously connected to the identity — which we 
call Lj 3 . We also pick up a 4 x 1 matrix "constant of 
integration" (C Q ) in obtaining (X a ) , so we are finally 
left with 



(5. 



For simplicity, we shall take this lowest-order piece to 
be the expected translation due to the position of ni2 at 
t = 0, viz., 



LcP — S a @, (Cq)q 



m 2 



(5.9) 



In fact, one can show that matching through third or- 
der requires that the spatial part of (C a ) be as given 
above (though the temporal part can still be freely spec- 
ified). Similarly, that matching requires Lj to differ 
from the identity only by a possible rotation about the 
y-axis, along with possible spatial and temporal reflec- 
tions. These are combined with a rotation that takes +y 
to — y if we have an odd number of reflections. We thus 
conjecture that matching at higher orders will further 
constrain this lowest-order coordinate transformation to 
be as given above. At the very least, the matching at 
fourth and fifth orders is independent of the remaining 
freedom. 



Similarly, noting that (bk) = 0, 



(9a 



P>2 



^(^ + 2 -6«-VWo} 



+ ^{3[(ri) -(S) Q ] 2 -[(r 1 ) ] 2 } 

+ ^{3[(n) ] 2 [(ri) • (S) ] - 5[(ri) • (b) f} 



A aj 3. 

(5.13) 



Here A Q/ 3 := diag(l, 1, 1, 1) is the "lowered 4-dimcnsional 

Kronecker delta." Also, we have (fik) = x k and, with 
our choice for the zeroth order coordinate transformation, 
( r i)o = £ ' equations for the coordinate transfor- 

mation at this order are thus 



[g a p) 2 = (^ 7 ) (M 2 (V)„ + ( A ^)i (Mo (V)r 

+2(A«r)„(M„(V) 2 . 

(5.14) 



C. First Order [0{^/rmJb)] 



since (h a b) 1 = 0, or, using our previous results, 



We have (g a p) 1 
matching gives us 



[h a p) l = 0. Thus, the first order 



(A«/3))i = °> 



(5.10) 



using the fact that (h a p) = r\ a p- As was given in 
Eq. (|5.6[) . the most general solution of this equation is 



(X a ) 1 = (Fpa)^ + {C a ) 1 



(5.11) 



where wc have written this in te rms of x a := x a — 
(m2/m)bx a because f/b := \J XkX k jb is one of our small 
parameters. Wc can do this without loss of generality, as 
it simply entails a different value for C a . 

This result differs from that given in Eq. (21) of Pa- 
per I; the latter suffers from some sign errors introduced 
during transcription. However, this does not affect that 
paper's final coordinate transformation, as the relevant 
constants were all taken to be zero. This is appropriate 
for Papers I and II, since they were using a corotating co- 
ordinate system: The boost encoded in our (F a p) 1 [seen 
in Eq. I|5.25[) ] would thus not be expected to appear in 
the coordinate transformation. 

D. Second Order [0(m 2 /b)] 

Proceeding to the next order, we have, recalling that 
(4)0 = 0, 



P>2 



2(Mi) b 



( gfc ')o ( X k ) (X 1 ) 
m 2 (R) b 2 1 > oK )o 



(5.12) 



(fla/9) a = [Kp) 2 + (F^), {F M \ + 2 . (5.15) 



Thus, at this order, the source of the differential equation 
[cf. Eq. (|5.4j) ] is given by 



2(5" Q ^) 2 = 2 (A {a/3) ) 2 = {g af }) 2 -{K f} ) 2 -{F a 1 ) 1 (F / g 7 ) 1 . 

(5.16) 

We now apply the integrability condition from 
Eq. (|5.5p and focus on the nonpolynomial piece of (S a p) 2 ; 
here this is the one that diverges [as (R) = (ri) = 
f — > 0]. It must satisfy the integrability condition in- 
dependently of the other pieces, by linear independence, 
and [considering, e.g., (X fci0 o) 2 ] g ivcs ( M i) = mi, as 
expected. The polynomial piece of the integrability con- 
ditions tells us that 



(£fci) = Ski - 3x k xi, (£ki P ) = 15x k xix p - 9<5( fci i p ), 

(5.17) 



using linear independence to read off the quadrupole and 
octupolc tidal fields separately. Solving for the coordi- 
nate transformation, we obtain 
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{X a ) 2 — 



1 - 



Aa/3i P + A ^ xl i a - I (F Q 7 ) 1 (^ 7 ) 1 ^ + (F/3q.) 2 i' 3 + (C„) 2 



(5.18) 



E. Third Order [0([m 2 /b} 3/2 )] 



At this order, the inner zone metric is 



to 2 V-R/i m 2 (i2) 

- COo (* fc ) (* ; 



?; 2 



b 3 



( g fc'p) i 

3& 3 



(5.19a) 



(**), = (*')„ (^)o + (*0o ( XP )o (^ s )o - (*')„ (* P )o (^)o • 



(5.19b) 



(^fcz) 3 = (^Oo) 3 $kl, 



and the near zone metric is 



m-2 V ? 'i / 1 & 



(5.19c) 



(5oo) 3 = ^ fr) -l[(^)o-Wi]+|{3[(n)o-Wi][(rl) -(6)o]-(n)i-(n)o}+^[(' 



&3' 



= 4?M 1 - TlV - iK^Oo ' Wo] + ^I^Oo ■ Wo! 2 - {(^i)o} 2 ] 



to [ (n) & 



2^[3{(ri) } 2 {(rl) • (&)„} - 5{(ri) • (6) } 3 ] 



(5.20b) 



(#fcz) 3 — (500)3 4z- 
I 



(5.20c) 



Here we have (b) 1 = ^Jmjm^tij, so 
— (7712/771) (6) 1 = — \J rri2 1 mty , and thus (ri) 1 • (&) 
0. Additionally, (6)j = \p^Jrn%{tjV)y and (&fe) 1 
y 7 m/rri2yk- We also have 

1 (A fc ) (A fc ) 1 



and 



1 

r i/ 1 



1 (n) • ( f 0i _ /W 



(n) [(n) ] 2 



(5.21) 



— ^3, (5.22) 



where f a := r a — (777,2 / l m)bx a . 

At this order, the equations for the coordinate trans- 



formation are 

2 (S a p) 3 = 2 (A( a ^)) 3 = (5^)3 - (^0/3)3 

- 2 (fc(a| 7 ) a ~ 2 (^l 7 )! (^1/3)7)2 . ( 5 -23) 

recalling that (h a p) 1 = and utilizing our lower-order 
results. To obtain (F a p) 1 , we look at the f~ 7 piece of 
(2fc;oo)3- Such a piece can only come from two spatial 
derivatives both acting on r ~ 3 in the f ~ 3 pieces of (Soo) 3 ; 
those pieces, in turn, only come from (1/R) 1 and (l/ri) 1 . 
Therefore, using Eqs. (|5.2ip and (|5.22[) . the intcgrability 
conditions require that we have 



[(F ak ) 1 x a + (C k ) 1 ] = -<r^yt, 

V TO 



from which we immediately see that 



m 2 * 



(5.24) 



(5.25) 
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using the antisymmetry of F a p. (Recall that to = — 1.) 

piece of (2^00)3 only comes 

0. 



By similar logic, the f 5 

from the (M 1 ) 1 piece of (^00)3 and gives us (Mi) 



The remaining nonpolynomial pieces cancel, so we have 
extracted all the information we can from the nonpoly- 
nomial part of the integrability condition. 

From the polynomial part of the integrability condi- 
tions, we first read off the octupole parts of the tidal 
fields, which arc 



6 / m 



— 9 / TTh 



and then the quadrupolc parts, which are 



— — / 7TZ 

(£ W ) 1 = 0, (B kl ) = -6J—x {k z iy (5.27) 



The third order coordinate transformation is thus 



m f yt 



x^x* 1 — Ax 2 



mjl b 2 ap I 2b 2 

T + ff- 2 )! + 5(-)' + 4 

oer V m Jo z \ m J 



(Ca)j 



1 

263 



xy(4x 2 - y 2 ~ z 2 )i a , 




r 



1 m 2 \ m 2 



yt a + 2 



1 - 



m 2 







2 771 / m 



(5.28) 



where the octupole part is the final term. 



Fourth and Fifth Orders [0([m 2 /b] 2 ) 
0([m 2 /bn 



and 



The matching at fourth and fifth orders proceeds in 
the same way as it did at lower orders, though the alge- 
braic complexity increases substantially. We shall thus 
give far fewer details of the calculations than we did be- 
fore, and mostly concern ourselves with pointing out the 
new features of the calculation that arise at these orders. 
The most prominent new feature, and the one responsi- 
ble for much of the algebraic complexity, is the presence 
of a nonpolynomial part in the coordinate transforma- 



tion. We know to expect this at fourth order because the 
transformation between Cook-Scheel and PN harmonic 
coordinates is nonpolynomial, and its lowest-order piece 
is 0(v 4 ) — see Sec. IV Gl However, there are various other 
nonpolynomial pieces present in the coordinate transfor- 
mation at fourth and fifth orders. We have to solve for 
these nonpolynomial parts of the coordinate transforma- 
tion by inspection (though we can still use Maple to 
obtain the polynomial part). It is reasonably easy to do 
so if one first breaks the source into pieces by multipolar 
order. 

The other subtleties involve the multipole expansion 
and are best illustrated by giving two examples: First 
looking at the fourth order piece of the inner zone metric 
component, we have 



(M,) 2 b I M) J1\ o [(Mx) f b 2 (£«) 2 / v * Wv in (*«)cwv*wv^ 



- 2 {x% {x% o + 2 (^H {x% {x% + §(^JWj M, (x\ (x% 



r 



The terms to notice are the ones involving (X fe ) 20 and £ki p /R- The first of these reflects the necessity 



22 



of avoiding "hidden octupole" pieces (i.e., pieces of 
octupolar order that arise when multiplying together 
pieces of lower multipolar order) when looking at cor- 
rections to lower-order terms: We only want to in- 
clude (X fc ) 2Q , the monopole (zeroth-order-in-x/6) piece 

of (X k ) r in (£ kl ) Q (X k ) {X l ) 2 Q /b 2 ; the dipole contri- 

bution to (^ fc ) 2 would give an octupole contribution to 



/ioo- Contrariwise, we do not want to leave out any terms 
of quadrupolar order, even if they arise from, e.g., oc- 
tupolar tidal fields. The £ki P /R contribution is such a 
term. 

These subtleties arise in a slightly different form in 
the equation for the fifth order piece of the coordinate 
transformation: 



2 (S afJ ) 5 = 2 (A [afi) ) s = (g a(} ) 5 ~ (h aP ) s - 2 (h {ah ) 4 (F^ - (i^) x (Ma (^A - 2(ftf a | 7 ) 3 {A W % 

~ 2(^) 3 0V) 7 ) 2 - 2 (F H i) 1 (h P rS ) 2 (A lP) s ) 2 - 2 (f (a| 'r) i (/^ _ (A w % fi (5.30) 

- 2(^ 7 ) 2 0V) 7 ) 3 ,< 3 - 2 (A Hl ) 2 (A W \^ - 2 (A (ah ) 2 (A m \ 2 - 2 (F (q , 7 ) 1 (A ]f} f) 4 . 



Here we have to avoid "hidden octupole" terms in many 
of the contractions: For instance, we only want the 
quadrupole-and-lower pieces of (-A( Q | 7 ) 2 (-A^) 7 ) , which 

are given by (A {ah ) 2 (A^f)^ + (i4(a| 7 ) a ,o (Atf) 3,2' 
However, it is also necessary to split up the inner zone 
metric's contributions into polynomial and nonpolyno- 
mial parts to keep from excluding quadrupole pieces as 
well [since the nonpolynomial parts of the second and 
third order pieces of the inner zone metric are all 0(b/R)]. 
In fact, this behavior means that we need to include the 
octupole part of the third order piece of coordinate trans- 
formation [in (/iSJL) 2 " y ) 3 <3 ]- (One would also need 
to include the octupole part of the second order piece of 
the coordinate transformation, but it vanishes.) 

The rest of the calculation proceeds as before, with 
the same general results: We obtain the next two orders' 
contributions to the matching parameters, with the non- 



polynomial pieces giving (F a p) 2 = and 

I / \ 3/2 / Tn 



+ 3 /^ +5 [W 
\ m / V m V m 2 



i[ a y/3], (5-31) 



along with (C fe ). = and (Mi). = for j E {2,3}. 
Similarly, the polynomial pieces give 



(£m) 
(B kl ) 



1 

2 ~ 2 



3x k xi - S k i H (4x fe x; - 5y k yi + hh) 

m 2 



5[^ 

m 2 



3/2 



7 HL-3 Tl 

y 1712 V 171 



(5.32) 



along with {Bki) 1 = {£ki) 3 = 0. The fourth order piece 
of the coordinate transformation is 



{X a ) i = -{At) i t ct + -^~ 

2 7712 



5_ _ 6- - - 

bf f b 



Tat ~l~ (>^a:)4 %a ~t~ (.Ay) ^ Da ~t~ (-^2)4 
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The fifth order piece is 
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Here t a &~i& is the 4-dimcnsional Levi-Civita symbol, with €0123 = 1- Additionally, for notational convenience, we have 
defined 
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where Mo (^4, B) and /C vanish in the limit 7771 — » (=> m — > 7772) 



G. Summary of matching results a relation between the mass parameters M and 777 1 , which 

we found to be equal [up to uncontrolled remainders of 
The final output of the matching is a set of expressions 0(v 4 )}; and the coordinate transformation necessary to 
for the tidal fields, which are given explicitly in Eqs. (|B1[) : 
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place the inner zone metric in the same coordinate system 
as the near zone metric to the order we matched. To 
obtain the full coordinate transformation, we start from 
Eq. (|5.ip and insert the various pieces we have read off 
or chosen. We shall take anything we were unable to 
fix by matching to be zero — here this will be Co (to all 
orders), along with (F a p)j and (Ck)j for j £ {4, 5}. With 
the results of our matching, this means that we have 
C a = (to all orders). We took the zeroth-order piece of 
the coordinate transformation to simply be the expected 
translation of the origin from the binary's center-of-mass 
to hole 1, so 

(X a ) Q = x a := x a - (m 2 /m)bx a . (5.38) 

With (C a )i = 0, we also have 

{X a ) 1 = (F f3a ) 1 S:' 3 , (5.39) 

where 



(5.40) 



Continuing onward, (X a )„ is given by Eq. (|5.18[) . where 
(F a p) 2 = and (C Q ) 2 = 0; {X a ) 3 comes from Eq. (pT28]) . 
where 



(F a p) 3 = M(0,5,3,l)t [a y m , 



(5.41) 



and (C a ) 3 = 0. [M(A, B, C, D) is defined in Eq. ([537]) .] 
Similarly, (X a ) 4 can be obtained from Eq. (|5.33p . and we 
take {F ol p) i = and (C Q ) 4 = 0. Finally, (X a ) 5 is given 
in Eq. (|5.35[) ; again, we set (F a p) 5 = and (C a ) 5 = 0. 

Looking back over this coordinate transformation, it 
is possible to gain some physical intuition about what 
it is accomplishing: The expected Lorentz boost due to 
the holes' orbital motion is present (through third or- 
der, which is the highest order at which we have fixed all 
the coordinate transformation, up to a possible temporal 
shift). We also have the lowest-order piece of the trans- 
formation between Cook-Scheel and harmonic PN coordi- 
nates for an unperturbed Schwarzschild black hole. [This 
is given by the first term in (At) 4 in Eq. (|5.34|) .] The 
remainder of the coordinate transformation is probably 
mostly concerned with effecting the transformation from 
locally inertial coordinates centered on the black hole to 
PN baryecntric coordinates. In addition, as we noted 
previously, the polynomial and nonpolynomial parts of 
the full coordinate transformation are related to the in- 
dividual holes in the expected manner: The nonpoly- 
nomial parts are associated with hole 1 and vanish in 
the limit mi —> 0. The polynomial pieces are associ- 
ated with hole 2 — indeed, everything except for the piece 
of the background Cook-Schccl-to-PN-harmonic transfor- 
mation vanishes in the limit to 2 -—> 0. 



4-metric g a p (in harmonic coordinates) in both the near 
and far zones. The resulting far zone metric is expressed 
in terms of derivatives of multipole moments obtained 
by integrating the "effective" stress-energy pseudotensor 
over the near zone. One also obtains nonlinear contribu- 
tions from integrating over the far zone (known as the 
outer integral in the DIRE approach), though only two 
of the resulting terms appear in the metric perturbation 
h al3 [defined in Eq. (2.2) of [6(|] to the order we are con- 
sidering. 

In this formalism, the metric perturbation in the far 
zone can be expressed in terms of the source multipolcs 
I® and J® via [Eqs. (5.12) in [U] 7 
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Here r := ||a?|| is the distance from the binary's center-of- 
mass to the field point and h k := x k /r is its associated 
unit vector. The (2"/r) 2 terms are the two contributions 
from the outer integrals mentioned previously We have 
included all the terms Pati and Will give, even though 
some of the ones in the purely temporal and spatial com- 
ponents only contribute terms that are of a higher order 
than we need here. We do this for completeness and also 
because we shall need these higher-order terms when we 
construct the extension to this data set (in Appendix [Cj. 

The source multipoles and J® (Q is a multi-index) 
are defined in Eqs. (4.5) of [HO] With these definitions, 
the mass monopole 2 is simply (a PN corrected version 
of) the total mass of the system, the dipole moment I k is 
the center of mass vector (so it vanishes in our coordinate 
system), and the current dipole J k is the total angular 
momentum. One can show that these three quantities 
are conserved up to radiative losses. 



VI. FAR ZONE METRIC 

7 We have corrected a sign error in Pati and Will's expression for 
h kl : They give coefficients of +2/3 and +8/3 for the second and 

The direct integration of the relaxed Einstein equations third terms . The correct signs can be obtained from Pati and 
(DIRE) approach [60,[61( can be used to compute the full Will's Eqs. (2.13) and (4.7b) in [g(|. 
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The source multipoles can be expanded in the PN ap- 
proximation to find 
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(6.2a) 
(6.2b) 

(6.2c) 

(6. 2d) 
(6.2e) 



The first two come from Will and Wiseman's 
Eqs. f4.16 ) |71j . and the remainder from Pati and Will's 
Eq. (Dl) [61|. (Even though we only need the lowest- 
order piece of X kl here, we include the 1PN correc- 
tions that Pati and Will give since we will need them 
in our construction of the higher-order extension in Ap- 
pendix [C]) Here the notation is mostly the same as for 
the near zone metric and was defined in Sec. IIV1 One 
new definition is x k := x\x\ (with similar definitions 
holding for different vectors, as well as larger collections 
of indices). 

Even though this structure is somewhat obscured in 
the expression Pati and Will give for it, the far zone PN 
expansion of the metric perturbation is clearly a bivari- 
ate expansion. There is a post-Minkowskian expansion in 
powers of G, given by radiative multipole moments, and 
a post-Newtonian expansion in v (or, more formally, 1/c) 
of these radiative multipole moments in terms of source 
moments computed in the near zone (see, e.g., [IH). Al- 
though formally the near and far zone expansions are 
independent, in practice there is a relation between their 
expansion parameters. This is given by the definition of 
the far zone, r > A = b/2v, from which one finds that 



since m/r comes with a factor of G due to the post- 
Minkowskian expansion while X\/r does not. This also 
means that we treat the (I/r) 2 oc (m/r) 2 term as 0(v 4 ), 
not 0(v 5 ). 



With this order counting in mind, we now determine 
the order of each source multipole moment. This leads 
to the following orders for the components of the metric 
perturbation in the far zone: h oa = 0(v 2 ), h ok = 0(v 4 ), 
and h kl = 0(v A ). Note that this is not the standard 
order counting of the metric perturbation in the near 
zone, since there we have b/r = 0(v°), which thus leads 
to h 00 = 0{v 2 ), h Qk = 0(v 3 ), and h kl = 0(v 4 ). The far 
zone order counting allows us to expand the full 4-metric 
in the far zone to obtain [from Eqs. (4.2) in [ic" 
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with remainders of 0(w 6 ). Here we have kept the 
— (l/2)h 00 term in the expression for <?ofc for formal con- 
sistency. It only gives terms of 0(v 6 ) with our order 
counting (which we neglect here), but would give terms 
of 0(v 5 ) if we had used the standard near zone order 
counting, where h ok = 0(v 3 ). This is the only 0(v 6 ) 
term in the expression for g fci so the uncontrolled re- 
mainder in that expression is thus actually 0(v 7 ). 
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(6.3) 



As explained in Sec. IV C of [60j, formal consistency re- 
quires that we treat each additional inverse power of r 
as raising the effective order by at least one factor of 
v. [This order counting is also equivalent to the order 
counting Alvi [4~fl ] uses in computing the far zone metric, 
though he keeps terms through 0(v 5 ) in all components.] 
For example, with this far zone order counting 
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We can now easily compute the full metric by carry- 
ing out all the differentiation in the expression for the 
metric perturbation, Eqs. (|6.1ap . In doing this, it is 
important to realize that the multipoles depend on re- 
tarded time, which must be carefully accounted for when 
differentiating. Eqs. (|6.5[) then give the full metric. In 
performing this calculation, we assume a (quasi)circular 
orbit, so the acceleration and the trajectories are paral- 
lel or antiparallel to each other and the velocity is per- 
pendicular to either of them to 0(v 5 ). We thus have 
xi-vi = 0(v ) = Si-Vi and a a = —uj 2 x*a + 0(v 5 ), where 
a A : ~ "A i s the acceleration of point particle A. [The lat- 
ter leads to expressions such as (v\ • ji) = — uj 2 v\ + 0(v 5 ), 
where j\ := a k A is the jerk of point particle A.] After 
much algebra, we finally obtain the full metric in the far 
zone: 



2G 



.9oo + 1 



9ki - S. 



kl 



2m\ mi f 2 m 2 , n f-, 2m | (x*i ■ n) , M x\ (x\ ■ n) jo 2 2^ 



- — + 2 fa • n) - — + 6- 
o r 



fa • n) - 4- + 



3 - 2r z uj' 



(6.6a) 



+ (1^2)+0(?; 5 ), 



5ofe = 



in i 
r 



4fa • n) + 4 fa • ri) 2 - 6 (xi • nf lu 2 + 4 V - 7 + 12 V - ' fa ■ n) - 2 -±- + 6 2 



+ ^iL (£1 ■ n) + 8 (xi • n) fa • n) - 2 ^ + 6 (xi ■ ^ l^x* + (1 ^ 2) + 0(v 6 ), 
r |_ r r J 



2mi mi 
On H 

r r 



2 m 2 



vi - — + 2 fa -hf + ^ +6 



to (x*i ■ h) 



+ -n kl + 4uf - 4w 2 xf I + (1 <-► 2) + 0(v & ), 
r 



r r 



(tTi • n) - 4 + (3 - 2^ 2 r 2 ) 



(6.6b) 



(6.6c) 



where everything is evaluated at the retarded time u := 
t — r. [This expression agrees with Alvi's result, given 
in his Eq. (2.17) [il|, though he also includes the 0(v 5 ) 
terms in the purely temporal and spatial components. 
We include these pieces — as well as even higher-order 
ones — when we construct the extension to the data in 
Sec. I VIII and have checked that we agree with Alvi about 
all the 0(v 5 ) terms.] 

A. Evolution of the binary's phase and separation 

Even though the effects of radiation reaction on the 
binary's orbital separation and phase arc formally small, 
only beginning at 0(v 5 ), they produce large corrections 
to the far zone metric when one is far away from the 
binary, since the retarded time at which one is evaluat- 
ing the binary's parameters becomes large. For instance, 
even as close as r — 50m, which is inside the outermost 
extraction radius (usually well inside) for all the simula- 
tions used in the Samurai project [12j, the phase differ- 
ence between a circular orbitfusing the 3PN expression 
for u> given in Eq. (190) of [64I] and 3.5PN inspiral (the 
computation of which is detailed below) is ~ 0.015 radi- 
ans for an equal-mass binary with a separation of 10m. 
(This phase difference should be compared with the av- 
eraged frequency domain phase accuracy required for pa- 
rameter estimation with Advanced LIGO, viz., 0.007 ra- 
dians, from p^ |. While such a comparison is not really 
warranted — see fl5| for some discussion — it gives a rough 
idea of the required accuracy.) See Sec. IV A in Kelly el 
al. [43[ for further discussion of the necessity of using for- 
mally higher-order PN results in obtaining the far zone 
metric. 

We thus use the most accurate (3.5PN) expression for 
the inspiral, as given by Blanchet [62| , who obtains it 
from an energy balance argument. The phase itself is 
given by Blanchet 's Eq. (234) [62|], where it is expressed 
in terms of the dimcnsionlcss time variable 0, defined in 



his Eq. (232) [62| . Since we are evaluating everything 
at the retarded time it, we have G = (r)/5m)(t c — u). 
Here the "coalescence time" t c is defined as the time at 
which the binary's frequency goes to infinity (or, equiv- 
alently, its separation goes to zero). One can calculate 
this in the PN approximation by using the energy bal- 
ance relation dE/dt = — C, where E and C are the bi- 
nary's energy and gravitational wave luminosity, respec- 
tively. These are given in terms of 7 := m/b though 
3.5PN in Blanchet's Eqs. (191) and (230) respec- 
tively. We can then compute t c by integrating dt/db = 
(dE/db)(dt/dE) = -{dE/db)/C from b = to b = b 
(where 60 is the binary's separation at t = 0). Here we 
expand the quotient as a power series (to 3.5PN) rather 
than using a Pade approximant (or performing any re- 
summation of the energy or luminosity) , as is sometimes 
done in the literature (see, e.g., [92|). 

With the 3.5PN expression for t c in hand, we can sim- 
ply substitute it into Blanchet's Eq. (234) [62| to obtain 
the phase as a function of u, making sure to expand to 
3.5PN order after substituting. We add a constant to 
the phase so that it is zero when t = 0, to be consis- 
tent with our choice of initial phase in the matching (i.e., 
so the holes initially lie on the x-axis, with our expres- 
sion for the orbit). We also take the freely specifiable 
gauge constant r and constant of integration 80 to be 
to and 1, respectively. (The dependence of u on u is 
then obtained by differentiating the phase with respect 
to time.) To obtain the retarded time dependence of 
6, we use b = m/7, along with the expressions for 7 in 
terms of x [Blanchet's Eq. (193) [HI] and x in terms of 
[Blanchet's Eq. (233) [62|]], expanding the quotient con- 
sistently to 3.5PN. (We take the appearance of b in a 
logarithm in the expression for 7 to consist only of 6's 
lowest-order dependence on 0, viz., 4TO0 1 / 4 .) N.B.: The 
final expansions of the expressions for (f> and b are impor- 
tant. If they are not performed, then one does not recover 
the expected values for b and uj at u — [viz., bo and u>o, 
respectively — here ujq is the binary's 3PN angular veloc- 
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ity for b = b obtained from Blanchet's Eq. (190) [62]]]. 
We do not display the resulting expressions, as they are 
quite lengthy, and best handled entirely within a com- 
puter algebra system. (We have carried out the calcu- 
lations in Maple and our scripts are available at [36l].) 
Our results can be seen graphically in Fig. |4j where we 
plot the past history of an equal-mass binary's separa- 
tion (starting from &o = 10m) along with the fractional 
deviations of its phase from lo$u. 




-300 , -200 
u/m 



FIG. 4: Fhe 3.5PN results for the past history of an equal- 
mass binary's separation, starting from 10m at u = 0, along 
with the fractional deviations of its phase from ujqu. 

Since we are effectively using higher-order equations 
of motion, due to the higher-order phasing relations, we 
would have liked to include higher-order terms in the 
relative-to-center-of-mass (relative-to-COM) variable re- 
lation as well, for consistency. However, the resulting 
expressions for the far zone metric components are alge- 
braically too complex for Maple to handle, so we had 
to forego including these terms. If one expands the con- 
tribution of the relative-to-center-of-mass relation to the 
(far or near zone) metric components consistently, then, 
to the order we have considered, all that one needs is the 
lowest-order Newtonian version of this relation, given in 
Eq. (|4.3[) . and this is all we have used in the far zone 
metric. Fortunately, unlike the secular radiation reaction 
effects in the binary's phase and separation considered 
above, the neglected PN corrections to the relative-to- 
COM relation are numerically small in addition to be- 
ing formally small: They first enter the metric at 0(v 8 ) 
(in the purely temporal component). Moreover, max- 
imizing over mass ratio (all the PN corrections to the 
relative-to-COM relation vanish for an equal-mass bi- 
nary), the magnitude of the largest of these is ~ 5 x 10~ 7 
at the inner edges of the near-to-far buffer zone for a 
separation of b = 10m; for comparison, the uncorrected 
(mi/r)vf + (1 <-> 2) term is ~ 10~ 4 in this situation. 

We were able to include these additional terms in the 
relative-to-center-of-mass relation in the near zone met- 
ric: For consistency, we include the same higher-order 
terms in the phasing and separation in the near zone 
metric as in the far zone metric, so it makes sense to 



attempt to include the higher-order relative-to-COM re- 
lation there, even though we were unable to do so in the 
far zone metric. In fact, these corrections contribute to 
the near zone metric at one order lower than to the far 
zone metric [viz., 0(v 7 ) in the spatiotemporal compo- 
nents] and are numerically quite a bit larger: Maximiz- 
ing over mass ratio, the magnitude of the largest of these 
corrections is ~ 2 x 10 -3 at the inner edge of the inner-to- 
near buffer zone for a separation of 10m; for comparison, 
the magnitude of the uncorrected — 4(mi/ri)i;* piece is 
~ 0.2 in this situation. 

Blanchet and Iyer give these corrections through 3PN 
in Eqs. (3.11)-(3.14) of [H]. We have specialized their 
result to a circular orbit by using the 2PN expression 
for to to express v in terms of m, b, and rj. Blanchet 
gives this relation specialized to a circular orbit through 
2.5PN in Eq. (187) of [62j|, so we shall just quote the 3PN 
contribution to x a here: 



?3PN 



-'/- 



(mi — m,2)m 2 



b 3 



7211 



22 

V 7T log 



b. 



In this expression, 



1260 ' 3 " \r 

'(6.7) 
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constant which, though a prion 
shall take to have the same value, viz., m. [This is equiv- 
alent to taking the related gauge constants and r' 2 
to both be m— see Eqs. (3.15) and (3.19) in Q.] The 
expression for the binary's separation vector that we sub- 
stitute into the resulting relative-to-center-of-mass rela- 
tion to obtain the trajectories of the point particles is 
b = b(xcos(j) + y sin0), where b and 4> are functions of t 
(in the near zone) or u (in the far zone). 



VII. INCLUDING HIGHER-ORDER TERMS 

Here we construct an extension to our data using var- 
ious readily available higher-order terms. This extension 
includes all the 0(v 5 ) terms in the near and far zones, 
but also includes even higher-order terms that do not 
improve the data's formal accuracy. [We have to change 
our far zone order counting somewhat to be able to claim 
that we have all the 0(v 6 ) pieces of the spatiotemporal 
components of the far zone metric; these components are 
necessary to obtain initial data valid through 0(v 5 ).} The 
general philosophy is that adding higher-order terms can 
often improve the quality of the data in practice, even if 
it does not improve their formal accuracy. As we have 
seen in the previous section, this is particularly true in 
the far zone, where the binary's phase evolution depends 
sensitively on the inclusion of quite high-order radiation 
reaction terms. 

There are also more specific reasons for including cer- 
tain of these terms: We would like for a putative evo- 
lution of our data to be able to be compared directly 
with Kelly et aZ.'s evolution of their data [57[ • (Such a 
comparison will give an indication of how much of the 
junk radiation is due to the failure of the initial data to 
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Zone Attribute _„ „ „ Versions 

04_No0ct 04 05 all 

r Time dependence Perturbative Perturbative Full Full 
Inner s 

I Fourth order octupoles No Yes Yes Yes 

^ r | Metric order 0(v 4 ) 0(v 4 ) 0(v 5 ) Full extended 

TABLE II: An overview of the contents of the various data sets we considered. 



include the correct tidal deformations.) Kelly et al. in- 
clude the 0(v 5 ) pieces of the spatial metric in the near 
zone [though not the matching 0(v 6 ) pieces of the ex- 
trinsic curvature]. The extension we have constructed 
includes these terms, as well. It also includes (as noted 
above) the 0(v 5 ) terms in the far zone, along with the 
0(v e ) terms in the extrinsic curvature in the near and 
far zones, so the extended data are valid through 0(v 5 ) 
in those zones. However, we have coded our data in such 
a way that one can easily produce a data set that only 
includes the pieces that Kelly et al. have, or some other 
subset of the pieces in the extension. 

We would have liked to have included the 0(v 5 ) terms 
in the inner zone as well, for completeness, but it is not 
possible to obtain initial data in that zone that is formally 
0(v 5 ) while still including the quadrupole pieces: We 
would need the 0(v 6 ) pieces of (the spatiotemporal com- 
ponents of) the inner zone metric. The quadrupole parts 
of these include hexadecapolc tidal fields, and knowledge 
of how those fields enter the inner zone metric would 
require nonlinear black hole perturbation theory. How- 
ever, we are able to calculate the polynomial parts of the 
fourth order octupole pieces — the results are presented 
in Appendix IB 21 — and include them in the extension. 
(These pieces include the 1PN correction to the electric 
octupole tidal field.) 

Moreover, we are also able to include the full time de- 
pendence of the tidal fields: See Appendix IB 21 It seems 
desirable to include these terms, since some of them are 
necessary for obtaining the pieces of the extrinsic curva- 
ture that describe how the tidal fields evolve: We already 
have the necessary, linear-in-t pieces for the lowest-order 
quadrupole fields, but not for their 1PN corrections, the 
explicit appearance of their time derivatives, or any of the 
octupole fields. As is discussed more fully in Appendix |D| 
the additional time dependence we include in the tidal 
fields is only the lowest-order dependence on t alone, 
not the additional space-dependence that comes from the 
tidal fields' dependence on the Eddington-Finkclstcin co- 
ordinate: The neglected dependence would generate con- 
tributions to the metric that combine with the unknown 
contributions from time derivatives of the tidal fields, 
while the t-dependent contributions are not entangled in 
this manner. We use t instead of T since the full time 
dependence of the tidal fields comes, in effect, from per- 
forming the matching at different (near zone coordinate) 
times t. We thus want these higher-order contributions 
to the tidal fields to depend on that time, instead of T 
(the inner zone time coordinate). 



We have not attempted to obtain the full time de- 
pendence of the coordinate transformation, since this is 
a rather more involved task than obtaining that of the 
tidal fields. Moreover, our rationale for including the 
tidal fields' full time dependence was to improve the evo- 
lution of the tidal perturbations: Including the full time 
dependence of the coordinate transformation would only 
improve the agreement of the inner and near zone met- 
rics in the buffer zone, while the largest effect of the tidal 
perturbations (e.g., in reducing the high-frequency junk 
radiation) presumably comes from closer to the horizons. 

The other higher-order pieces we have added are all in 
the purely temporal components of the near and far zone 
metrics. We have included these since Blanchet, Faye, 
and Ponsot [58[ give an explicit expression for the purely 
temporal component of the near zone metric through 
0(v 7 ) and it is easy enough to calculate the matching 
terms in the far zone. The specifics of where we obtained 
all the extra terms are given in Appendix [Cj along with 
the accompanying caveats and new order counting in the 
far zone. We refer to these "state-of-the-art" versions of 
the near and far zone metrics as the full extended ver- 
sions. 

We have considered the effects of these additional 
terms on how well the metrics stitch together by putting 
together four versions of the data, as shown in Table [H] 
While the versions in the table are given the labels we use 
for the corresponding Maple scripts and C code (avail- 
able at [36[ ) , we shall usually refer to them by the order 
of their near and far zone metrics [i.e., as 0(v 4 ), 0(v 5 ), 
and full extended]. This will not cause confusion, since 
we will almost always consider just the right-most three. 
While the 04_No0ct version is important, since it con- 
tains only the pieces that can be included consistently in 
the matching (except for the higher-order terms neces- 
sary to obtain the phasing in the far zone accurately and 
the analogous terms in the near zone), we only consider 
it in Fig. 1211 Thus, when we refer to the 0(v A ) version 
without any qualifiers, we mean 04, the version includ- 
ing the fourth order octupole pieces. See Appendix [D] 
for the specific details of how all these metrics are cal- 
culated. (N.B.: The near zone metrics in all of these 
versions contain the background resummation detailed 
in Sec lVIIIAH 

In the next section, we compare the extended data with 
their lower-order counterparts: The volume elements of 
the 0(v 4 ), 0(v 5 ), and full extended sets of data are com- 
pared in the near zone in Figs. [5] and [5] and in the far 
zone in Figs. [7] and [SJ Their constraint violations are 
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compared in Figs. |T7l [TH and [T9j Wc compare the con- 
straint violations of the different versions of the inner 
zone metric in Figs. 1201 and I2T1 In general, the additional 
terms reduce the constraint violations and improve the 
overall matching, or at least do not affect them adversely. 
The only exception to this is the addition of the full time 
dependence in the inner zone, which increases the con- 
straint violations in the inner zone. However, it does not 
increase them by much, and they were originally quite 
small: Our philosophy is thus to include these terms, 
which wc think will improve the data's evolution, even 
at the cost of slightly higher constraint violations. 



VIII. NUMERICAL CONSIDERATIONS 

To get an idea of how well the matching is working 
numerically, we plot the volume elements of the various 
4- metrics in the vicinity of their buffer zones. (In all these 
plots, we use the full extended version of the data dis- 
cussed in the previous section, unless otherwise noted.) 
Of course, we first transform the inner zone metric using 
the coordinate transformation obtained from the alge- 
braic matching above. To obtain a simple, easily inter- 
preted plot, we choose the test system b = 10m, mi = m,2 
(the mass is in arbitrary units) and restrict our attention 
to the t = timeslice and the x-axis (i.e., we consider 
the spatial slice along the separation between the holes), 
concentrating on the portion near hole 1. (We expect this 
slice to provide the most stringent test of the matching, 
since it contains the portions of the buffer zones where 
the field is strongest and changes most rapidly.) We shall 
use this setup (or slight modifications thereof) for all of 
our later examples. The matching in the inner-to-near 
transition is shown in Figs. [5] and [5] (along with a compar- 
ison of different near zone volume elements). The first of 
these shows the volume elements themselves [along with 
the differences between the volume elements of the 0(v 4 ) 
and 0(v 5 ) versions of the near zone metric] and the sec- 
ond shows the differences between the volume elements 
of the various near zone metrics (along with the merged 
metric) and the inner zone metric. The volume elements 
in the near-to-far transition are shown in Fig. [7] (along 
with the volume element of the merged metric). 

The inner zone metric we display here contains the 
fourth-order octupole pieces discussed in Appendix IB 21 
To avoid clutter, we did not plot the "plain" version with- 
out these additional pieces: It agrees very closely with the 
version we have plotted near and inside the horizon, but 
bends away from the near zone metric further away from 
the hole. The differences between the versions of the 
inner zone with full and perturbative time dependence 
would not appear in this plot, since we are looking at the 
t = timeslice. We display the differences between the 
0(v 4 ), 0(v 5 ), and full extended far zone metrics in Fig.0 
considering a binary with mass ratio 3 : 1 so that these 
differences are more pronounced. [In this plot, we have 
concentrated on the portion of the x-axis to the right of 
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FIG. 5: The volume elements of the transformed inner zone 
metric (dotted line) along with those of the 0(v 4 ) and ex- 
tended versions of the near zone metric. These are plotted 
both with [0(v 4 ): dot-dashed line; extended: solid line] and 
without [0(v 4 ): dashed line; extended: dot-dot-dashed line] 
background resummation. We have also displayed the differ- 
ences between the volume elements of the 0(v 4 ) and 0(v 5 ) 
resummed metrics in the inset. We calculated all of these vol- 
ume elements for our equal-mass test system and have plotted 
them along the x-axis near hole 1. (The associated point par- 
ticle is at x/m = 5.) We have indicated the rough locations 
of the intersection of the x-axis with hole l's inner-to-near 
buffer zone (BZ) using double-headed arrows. 
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FIG. 6: The differences between the volume elements of var- 
ious metrics with that of the full extended inner zone metric. 
The various metrics considered are the near zone metrics plot- 
ted in Fig. [5] (using the same colors and line styles as in that 
figure) as well as the merged metric constructed in Sec. lVIIIBI 
We have indicated the rough locations of the intersection of 
the a;-axis with hole l's inner-to-near buffer zone (BZ) using 
double-headed arrows. 



the more massive hole, for clarity: The relations between 
the various metrics arc qualitatively the same to the left 
of the less massive hole, except for differences in relative 
amplitude.] 

In the near zone metric, the extension adds certain 
terms that become large when evaluated close to the holes 
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FIG. 7: The scaled volume elements of the (resummed) near 
zone [both 0(v 4 ) and extended], far zone, and merged met- 
rics. These are calculated for our equal-mass test system and 
plotted along the x-axis in the portion of the near-to-far tran- 
sition region that lies to the right of hole 1. To obtain a less 
crowded plot, we have scaled these by the contribution to the 
volume element of the lowest-or der nontrivia l pieces of the far 
zone metric, viz., (1 + 2m/r)\Jl — Am 2 /r 2 . The buffer zone 
is roughly the region in which we transition, which, in turn, 
is roughly the region plotted. 




x/m 

FIG. 8: The scaled volume element of the far zone metric 
for a binary with a separation of 10m and a mass ratio of 
3 : 1 along the x-axis to the right of hole 1 (the more mas- 
sive hole). We display this for the extended data (the darkest 
curve — black in color) as well as the 0(v 5 ) and 0(v 4 ) ver- 
sions of the data (the medium and lightest curves — blue and 
red in color, respectively) to illustrate their differences. (We 
considered an unequal mass binary to make these differences 
more pronounced.) As before, we have scaled these volume 
elements by the contribution from the lowest-order nontrivial 
piece of the metric, viz., (1 + 2m/r)^Jl — 4m 2 /r 2 . 

(e.g., Am\m2/r\ in goo, which equals 1/2 at r\ = m for 
our test system) and others that grow rapidly as one 
moves away from the holes (e.g., — m^m^rf /46 6 , also in 
goo, which is about —0.03 at n = 100m for our test 
system). The terms that become large near the holes 
cause the t = timeslice of the near zone metric to no 
longer be spacelike in a region that extends outside the 



horizon. This can be seen in the rapid decrease of the 
volume element around the hole in Fig. [5j (Even though 
we do not display this in the plot, the volume element in 
fact decreases to zero.) However, for this separation, the 
timeslice is still spacelike in the buffer zone, so its bad 
behavior closer to the holes does not cause any problems 
in the merged metric. 

N.B.: The unperturbed horizon is mi[l — vnijh + 
0(m^/b 2 )] away from the point particle associated with 
hole 1 in the new coordinates. For the test system, the 
correction is small, so hole l's unperturbed horizon inter- 
sects the x-axis at ~ 4.5m and ~ 5.5m in the new coordi- 
nates. There are also corrections due to the tidal distor- 
tion, but for the test system these are even smaller than 
those due to the new coordinates: Taylor and Poisson [73[ 
study the effect of the tidal distortion on the horizon (in 
their Sec. VIII) and find the expected quadrupolar de- 
formation at lowest order. 

In the inner-to-near transition, the volume elements 
generally behave in the expected manner: They differ 
when considered either too near or too far from the hole 
and approach each other in the buffer zone. Things look 
a bit more unusual in the near-to-far transition, since the 
two volume elements agree better in the region between 
20m to 30m than they do further away from the holes. 
However, this should not be surprising: The reduced 
wavelength of the gravitational radiation is ~ 16m for 
the equal-mass test system, and the difference between 
the near zone metric's perturbative treatment of retar- 
dation and the far zone metric's full treatment should 
become quite apparent beyond that radial distance. In 
fact, the oscillations we see in the far zone metric's vol- 
ume element in Fig. [5] are due to the far zone metric's 
dependence on retarded time. (Even though the met- 
ric contains gravitational radiation, this is not the only 
source of these oscillations. However, it does contribute 
to them, as expected.) This plot also illustrates the phase 
differences between the various versions of data. Addi- 
tionally, we expect the volume elements of the near and 
far zone metrics to differ dramatically close to the binary. 
This is indeed the case (particularly for r < 10m), though 
we did not show this in the plot, preferring instead to con- 
centrate on displaying the details of the volume elements' 
behavior in the transition region. 



A. Background resummation 

In the inner-to-near transition, the volume elements 
do not agree as closely as we might like. In fact, the 
agreement is worse near the hole with the higher-order 
extended version of the near zone metric than it is with 
the original 0(v ) version (though the agreement fur- 
ther away from the hole and outside the orbit is slightly 
better). However, even the agreement with the original 
version is not much better than that found in Paper II 
(see its Fig. 2), even though we have matched to higher 
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order. 8 Part of the resolution of this apparent problem 
is that we have not yet used one of the "tricks" from 
Paper II, namely background resummation. 

The idea of background resummation is to add higher- 
order terms to the near zone metric in order to im- 
prove its strong field behavior. If one considers the limit 
7712 — > (=>■ v\ — » 0) of the near zone metric, then one 
finds that it reduces to the far field asymptotic expansion 
of the unboosted Schwarzschild metric (with mass param- 
eter mi) in (PN) harmonic coordinates. This expansion, 
however, lacks the causal structure of the Schwarzschild 
metric — it has no horizon. 

One method to restore this causal structure is to re- 
sum the PN metric. This resummation consists of adding 
an infinite number of higher-order terms, such that the 
metric reduces identically to the full (not the expanded) 
unboosted Schwarzschild metric in (PN) harmonic co- 
ordinates in the limit ni2 — > 0. A priori, there is no 
reason to suspect that adding such higher-order terms 
would increase the accuracy of the PN metric. A poste- 
riori, however, it is usually the case that such rcsummcd 
metrics are indeed closer to the exact solution, as was 
seen in Paper II. Moreover, we should stress that these 
higher-order terms are not arbitrary, but guided by the 
physical requirement of restoring the apparent horizons. 
Such physically informed rcsummations have met with 
great success in general relativity, notably in the effec- 
tive one-body formalism (see, e.g., [HI). 

It would also be possible to resum the far zone met- 
ric: One can calculate its Newtonian part without mak- 
ing a multipolc decomposition by proceeding in the same 
manner as in the calculation of the Lienard-Wicchert 4- 
potcntial in electrodynamics (as given in, e.g., §14.1 in 
Jackson (9o|). However, we have not pursued this line 
of investigation further. Simone et al. [9t| performed a 
related resummation of the Newtonian pieces of the lu- 
minosity of an extreme mass-ratio binary, though they 
did this by first calculating the multipole expansion and 
then rcsumming directly, while the resummation of the 
metric we are suggesting here would involve computing 
the integral directly, with no multipole expansion of the 
Newtonian part. 

With these points in mind, let us now describe the pro- 
cedure in detail, exemplifying it using the purely tempo- 
ral component of the metric. In PN harmonic coordi- 
nates, this component of the Schwarzschild metric takes 
the form —(R — M)/(R + M). Here R is the harmonic 
radial coordinate and M is the hole's mass. (As our use 
of "i?" indicates, this is the same radial coordinate as in 
Cook-Scheel harmonic coordinates — see Appendix [Al for 



The figures are not directly comparable, since the one in Pa- 
per II plots the xx components of the metrics, not their volume 
elements. However, the plot of just the xx components of this 
paper's inner and near zone metrics displays the same behav- 
ior as that of the volume elements, including roughly the same 
numerical values for the difference between the metrics. 



a comparison of the two coordinate systems.) We expect 
to have mi — > M and r± — ► R as m<i —* 0, which suggests 
that the purely temporal component of the PN metric 
should approach — (ri — m\)j(r\ + mi) as — ► 0. The 
far field asymptotic expansion of this metric component 
(i.e., for n 3> mi) is given by 
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This identically reproduces all the terms in the PN near 
zone metric in the limit mi — > 0. We then resum the PN 
near zone <?oo by taking 
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where is the version of this component without re- 
summation, given in Eq. (|4.1a| . 

A similar procedure can be applied to the spatial sector 
of the metric. Carrying this out, we obtain 
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where := xf jr\ and g£] d is given by Eq. (|4.1cp . One 
can check that gki reduces to g%\ d identically as m 2 — > 0. 

We have used the Schwarzschild metric in PN har- 
monic coordinates to resum the PN metric here, since 
it is this background that the PN metric approaches in 
the ni2 — * limit. We thus cannot resum the spatiotem- 
poral components of the metric, since they already match 
the background in this limit. If we had first transformed 
the PN metric to Cook-Scheel coordinates, then we would 
have been able to resum the background so that it exactly 
matched that of the inner zone metric. This would have 
guaranteed a better agreement, and would have probably 
also given a merged metric with smaller constraint viola- 
tions, since there would be no coordinate singularity at 
the horizons in the resummed near zone metric. 

However, if we had chosen this route, we would have 
had to pick some region surrounding the buffer zone 
for each hole in which to perform this transformation. 
This would have introduced further complications that 
we thought it best to avoid in this implementation, even 
at the possible cost of somewhat poorer matching. More- 
over, the background resummation procedure that we 
have implemented has indeed improved the matching, as 
can be seen in Fig. [5j The improvement is particularly 
striking for the extended version of the near zone metric, 
where the region in which the t = slice is no longer 
spacelike moves closer to the horizon and the graph of 
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its volume element now crosses that of the near zone in 
the buffer zone outside of the orbit. But resummation 
also improves the matching of the original data: The 
resummed version of the 0(w 4 ) near zone metric agrees 
more closely with the inner zone metric than does the 
unresummcd version. 



B. Transition functions 

We now turn to the process of stitching the inner and 
near zone metrics together numerically. It is necessary to 
interpolate between the various metrics in transition re- 
gions (located inside the buffer zones) in order to stitch 
the metrics together with no discontinuities. A simple 
way to do this is to use a weighted average of the metrics, 
where the precise way this average is carried out is de- 
termined by a C°° transition function Fab '■ — > [0, 1]. 
This function should have the property that Fab (x) = 
if x e C A n AB and F A b{x) = 1 if x G C B D O ab - 
Here Oab * s the complement of the buffer zone between 
zones Ca and Cg. (These conditions guarantee that the 
transition takes place completely inside the buffer re- 
gion.) If we just consider two metrics, g^j and g^ , 
for simplicity, then the resulting merged metric is given 

by [1 - FAB]g a p + FAB9ap- ( Wc have suppressed the 
position dependence of the transition functions and met- 
rics, for notational convenience.) 

In principle, these conditions arc all one would impose 
on possible transition functions. (One might also want 
to stipulate that Fab be "increasing as one moves from 



where g^p denotes the metric that lives in zone A (see 
Fig. [5] for the numbering system) and 

/far(r) :=/(r,A/5,A,l,2.5), 
fneax(x) ■= /(x,2.2m 2 - m-^b/m, b - 2.2m, 1,2.5), 

/inneMM := f(r A , 0.256rX S.Uim'b 5 ) 1 / 7 , 0.2, b/m). 

(8.6) 

[N.B.: We shall refer to the value of a parameter ap- 
pearing in a particular transition function by trans- 



Ca to (7b, " where this would have to be interpreted in 
some appropriate sense.) One could then contemplate 
minimizing some appropriate norm of the constraint vi- 
olations of the resulting merged metric (outside, say, the 
apparent horizons) over all possible transition functions 
satisfying these requirements (97j . Our approach will be 
significantly less ambitious, leaving a systematic study of 
transition functions to later work, probably that accom- 
panying an evolution of the data. 

However, as Yunes [1^] realized, one can exclude 
from consideration a priori any transition functions that 
would induce larger (formal) constraint violations in the 
data than those due to the uncontrolled remainders of 
the individual metrics. The so-called Frankenstein the- 
orems enunciated in [35| provide sufficient conditions on 
the transition functions to keep this from occuring: The 
first Frankenstein theorem tells us that the first and sec- 
ond derivatives of Fab must be 0(ca, c_b) and O(e A ,e B ), 
respectively, if the merged metric is to satisfy the Ein- 
stein equations to the same order as the individual met- 
rics. (Here ca represents the small parameters associated 
with zone A.) We have constructed our transition func- 
tions to respect the conditions of this theorem — we check 
this explicitly below. 

We shall take our transition region and functions to 
be spherically symmetric, even though the optimal ones 
would surely be distorted — neither the binary nor the 
holes are spherically symmetric. Moreover, we shall only 
consider one particular form for the transition functions, 
viz., the same form used in Papers I and II: 
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ferring that transition function's subscript to the pa- 
rameter's name, e.g., Winner, A '■= 3.17(m 2 6 5 ) 1 / 7 .] Here 
m is the Newtonian wavelength of the bi- 
nary's gravitational radiation. We have also used the 
"transition radius" rT, given by taking the leading orders 
of the uncontrolled remainders of the approximations in 
the inner and near zones to be comparable. This gives 
{m/b){r T A /b) A = (m A /rl) 3 , and thus rj = (m A b 5 /m)^ 7 . 
With these choices for the parameters, the transition 
functions satisfy the hypotheses of the first Frankenstein 
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where x(r, ro,w) := tan[7r(r — ro)/2w]. The full merged metric is thus 

a/3 = {!- /far(? , )}{/noar(a;)[/inncr,l(?'l)gi 3 / ] + {1 - /inner, 1 (^1 ) }s , aj] 



g a p = 1.1- JiB,r{r)}^J neaI {X)lJ ianeI A{ri)g^ - /inner, l^lJI^J 

+ [1 - fneaxix)] [/inner,2 fa)^ + {1 - /inner,2 fa)}^] } + /far(r)s£ 

I 



(4) 
0' 



33 



theorem [35|, given above: For instance, the nth spatial 
derivatives of / noar and /f ar scale as iu~ e " r oc b~ n oc v 2n 
and = 1/A" oc v n , respectively. Matters are a 

bit more complicated for /inner, a: Its spatial derivatives 
decrease rapidly as v — > because they involve sech 
with an argument that goes to infinity as v — ► (since 
•Sinner. A oc v~ 2 ) and sech is rapidly decreasing at infinity. 

We determined the parameters we use for the transi- 
tion functions by experimenting with different choices: 
We found that the values given in Eq. (|8.6j) produced the 
smallest overall constraint violations along the x-axis for 
our equal-mass test system (with b = 10m) of all the 
choices we tried. Except for the near-to-far zone tran- 
sition, which is completely new, these choices are very 
similar to those in Paper II. (In fact, despite the way 
it is written, our / noar agrees exactly with its analogue 
from Paper II, the function G.) The most important dif- 
ference is the scaling of the transition width, Winner, A- 
Both Papers I and II took Winner, A to scale with rj, so 
it depends on the system's mass ratio; in particular, it 
goes to zero as uia — > 0. (Note that our /i, m or,A corre- 
sponds to Fa in Papers I and II.) We have found that 
this choice for Winner, A results in large transition-induced 
constraint violations near the smaller black hole for un- 
equal mass ratios. These constraint violations appear to 
increase without bound as the mass of the smaller hole 
goes to zero. This is as one would expect, since the gra- 
dient of the transition function blows up as the transition 
width goes to zero, leadin g to an unbounded increase in 
constraint violations — see [35| and the above discussion 
of the Frankenstein theorems. 

We can obtain a better-behaved transition function by 
freezing the dependence of Winner, A on the mass ratio, as 
we have done here. (We also have a slightly different co- 
efficient of r\ than in Paper II, even for an equal-mass 
system.) The other differences between our /inner, A and 
Paper IPs Fa are a slightly different coefficient of the 
transition radius in ro,inner,Aj m addition to a new tran- 
sition radius that reflects the higher-order matching that 
we have performed here. (The transition radius, r\ : was 
slightly misleadingly called the matching radius and re- 
ferred to as in Paper II.) 

The choices for the transition functions' parameters de- 
termine effective boundaries for the various zones: These 
are given in Table HTT1 for our equal-mass test binary. This 
table gives both the formal boundaries (i.e., the numer- 
ical values of the boundaries given in Sec. HJ) and the 
effective boundaries (i.e., the boundaries determined by 
our choices of parameters for the transition functions). 
What we call the complete effective boundaries are deter- 
mined by the entire region in which we use a given metric, 
even if the coefficient of the metric (due to the transi- 
tion functions) is very small in a portion of the region. 
What we refer to as the practical effective boundaries 
are cut off when the coefficient of the metric becomes 
smaller in magnitude than 10 -4 [i.e., much smaller than 
the magnitude of the uncontrolled remainders, which are 
~ (b/m)~ 5 / 2 ~3 x 10~ 3 ]. 



Zone 


Formal boundaries 


Effective boundaries 
Complete Practical 


Inner 


r A < 10m 


r A < 17 'Am r A < 11.2m 


Near • 


r A S> 0.5m, 


r A > 0.985m, r A > 1.27m, 


r < 15.8m 


r < 119m r < 109m 


Far 


r > 15.8m 


r > 19.9m r > 30.4m 



TABLE III: The zone boundaries for our equal-mass test bi- 
nary (b = 10m). 



The effective inner zone boundaries given in Table IIIII 
are not quite correct: Since even the practical effective 
boundaries of the inner zones are greater than half the 
distance between the holes for b = 10m (as they are for 
an equal-mass binary, with our choice of transition func- 
tions, for b < 165to), one needs to introduce a third tran- 
sition function, here called /near, to effect the transition 
between the holes. (See Sec. VI B of Paper I for further 
discussion; this function is referred to as G in Papers I 
and II.) With our choice of parameters, / noa r cuts off the 
complete effective inner zone when the ^-coordinate is 
closer than 1.1m to (the x-coordinate of the point parti- 
cle associated with) the other hole. The practical effec- 
tive inner zone is cut off when the x-coordinate is 1.93m 
away from the other hole.) 

In /far, one might be concerned that even the practi- 
cal effective transition region extends well outside of the 
standard outer boundary of the near zone, viz., r ~ X. 
Indeed, it is quite possible that blending in the near zone 
metric in a region where its perturbative treatment of 
retardation is not warranted will introduce significant 
phase errors in the binary's outgoing wave train. How- 
ever, our choice of transition region is justified (at least 
for this preliminary construction of transition functions) 
by the (relatively) large constraint violations of the far 
zone metric at the inner edge of our transition region, as 
is shown in Fig. 1131 



Merged 

Near zone 

Far zone 




FIG. 9: The norm of the shift (\/f3 k f3k, with indices lowered 
by the metric under consideration) of the near zone, far zone, 
and merged metrics to the right of hole 1 along the axis be- 
tween the holes. As usual, these are all calculated for our 
equal-mass test system. 
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FIG. 10: The trace of the extrinsic curvature of the inner 
zone, near zone, and merged metrics near hole 1 along the 
axis between the holes. As always, these are all calculated for 
our equal-mass test system. 
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FIG. 11: The trace of the extrinsic curvature of the near zone, 
far zone, and merged metrics to the right of hole 1 along the 
axis between the holes. As usual, these are all calculated for 
our equal-mass test system. 



To illustrate the matching, we have plotted the vol- 
ume elements of each zone's metric along with that of 
the merged metric, all for our standard test system. (Re- 
call that we use the full extended version of the data in all 
plots unless otherwise noted.) We did this for the inner- 
to-near zone transition in Fig. [6] plotting the difference 
between the merged and inner zone metrics' volume ele- 
ments, for clarity. Similarly, the plot for the ncar-to-far 
transition was already given in Fig. [7] We also need to 
consider the matching of the shift ((3 k ) separately, since 
the volume element of the 4-metric equals the lapse times 
the volume element of the 3-metric [see, e.g., Eq. (E.2.25) 
in Wald [HI]- We display its norm (y/ j3 k [3k, with indices 
lowered by the metric under consideration) in the near- 
to-far transition in Fig. [9] The analogous plot in the 
inner-to-near transition looks similar enough to that of 
the volume element (Fig. O that we do not include it. 
For the same reason, we do not show the matching of 



the norm of the extrinsic curvature (yj K 1 pKi p , with in- 
dices raised by the metric under consideration) in cither 
transition region. The trace of the extrinsic curvature, 
K, behaves differently enough that we plot its match- 
ing in Figs. [10] and [XT] Moreover, the trace of merged 
metric's extrinsic curvature is an intrinsically interesting 
quantity: It tells us how much our data's slicing differs 
from maximal slicing (K = 0). 



C. Constraint violations 

The constraint violations provide a much more sensi- 
tive check on how well the transition functions are work- 
ing than the previous subsection's plots, in addition to 
giving a measure of the accuracy of the entire initial 
data construction. We compare the constraint violations 
(computed using BAM [98|, [99[ ) of the individual metrics 
to those of the merged metric for our equal-mass test 
system along the x-axis in Figs. [12] and [13] (We have 
checked that the constraint violations behave roughly 
similarly — and are not significantly worse — in the y- and 
z-directions.) For the norm of the momentum constraint, 
Mfc, we have chosen y/ M k Mk, with indices raise d by the 
(merged) metric. [See, e.g., Eqs. (14)-(15) in [2l[ for 
expressions for the constraint equations.] N.B.: As one 
can see by looking at the (norm of the) momentum con- 
straint, the plot of the constraint violations in the inner- 
to-near transition (Fig. [12]) does not cover the entire tran- 
sition region. However, the portion we do not show is not 
particularly interesting: The momentum constraint vio- 
lations of the near zone metric continue to decrease, and 
those of the merged metric rapidly become indistinguish- 
able from the near zone metric's momentum constraint 
violations.) 

The largest constraint violations in the merged metric 
outside of the horizons occur in the inner-to-near tran- 
sition regions. These are likely due in part to the near 
zone metric's large constraint violations near the horizons 
due to the coordinate singularity at the horizons. (Re- 
call that the horizons are approximately .5m away from 
the positions of the point particles in our test system — 
the effects of the new coordinates and tidal distortion are 
small.) However, since the only knowledge we have about 
the size of the uncontrolled remainders is their scaling 
with v (and thus b), the magnitude of the constraint vi- 
olations does not really provide us with a check that the 
uncontrolled remainders are in fact of the advertised or- 
ders. We can obtain such a check by considering how the 
constraint violations vary as the binary's separation (b) 
increases. As expected, they decrease at least as rapidly 
as &~ 5 / 2 (the expected scaling of the largest uncontrolled 
remainders) as the binary's separation increases. This is 
illustrated in Figs. [T4l and 1 151 in the inner and near zones, 
and in Fig. [THlin the far zone. 

For a separation of 10m, the relatively large constraint 
violations of the merged metric compared to those of the 
individual metrics are an indication that this separation 
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FIG. 12: The Hamiltonian constraint violations and norm 
of the momentum constraint violations of the merged, inner 
zone, and near zone metrics. These are plotted along the 
a;-axis around hole 1 for our standard equal-mass test system. 
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FIG. 13: The Hamiltonian constraint violations and norm of 
the momentum constraint violations of the merged, near zone, 
and far zone metrics. These are plotted along the x-axis in 
a region including the portion of the near-to-far transition 
region to the right of hole 1. As usual, this is done for our 
standard equal-mass binary test binary. 



FIG. 14: The Hamiltonian constraint violations around hole 
1 for an equal-mass binary with separations of 10m, 15m, and 
20m. For ease of comparison, we have scaled the a>axis by 6 
so that (the point particle associated with) hole 1 is always 
at the same position. In the inset, we zoom in to show how 
the inner zone constraint violations vary with b. 




FIG. 15: The norm of the momentum constraint violations 
around hole 1 for an equal-mass binary with separations of 
10m, 15m, and 20m. For ease of comparison, we have scaled 
the s-axis by b so that (the point particle associated with) 
hole 1 is always at the same position. 



is close to the minimum for which the hypotheses under- 
lying the data's construction are valid: For instance, as 
seen above, for this separation (and an equal-mass bi- 
nary), the two inner-to-near transition regions overlap 
between the holes, meaning that much of the inner-to- 
near transition is effected by / nea r — see Table [Hi] and the 
surrounding discussion. Moreover, as we have seen pre- 
viously, the t = slice of the near zone metric is not 
even spacclike at some points outside the holes' horizons 
for b = 10m. However, this does not adversely affect the 
merged metric with our choices of transition regions. If 
one tries closer separations, things are significantly worse. 
For instance, for a separation of 6m, the maximum con- 
straint violations are larger than those for 10m by a factor 
of 10 or more. Moreover, with our choices for the tran- 
sition regions, the merged metric contains some of the 



portions of the near zone metric where its £ = slice is 
not spacelike. 

In the near-to-far transition (as shown in Fig. I13[) . the 
relatively large constraint violations associated with the 
far zone metric are not unexpected at the distances from 
the binary at which we have made our transition: We 
have uncontrolled remainders of, e.g., 0{mb A /r 5 ) in the 
multipolar expansion of the far zone metric, so we ex- 
pect it to have (dimcnsionless) constraint violations of 
0(m 3 6 3 /r 7 ). For instance, at r = 30m, m 3 b 3 /r 7 ~ 
5 x 10 -7 , reproducing the order of magnitude of the con- 
straint violations at that point. We would have thus 
needed to transition somewhat farther from the binary 
if we wanted the near and far zone metrics' constraint 
violations be comparable in the transition region. How- 
ever, we almost surely do not want to do this, since the 
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FIG. 16: The Hamiltonian constraint violations and norm ol 
the momentum constraint violations along the a>axis in the 
portion of the near-to-far transition and far zone proper that 
lies to the right of hole 1. These are computed for an equal- 
mass binary with separations of 10m, 15m, and 20m. For ease 
of comparison, we have scaled the x-axis by A = n\Jb A /m, the 
Newtonian wavelength of the binary's gravitational radiation. 
In the two insets, we zoom in to better illustrate the behavior 
of the Hamiltonian constraint violation in two situations: In 
the lower inset, we consider the b = 20m data in the transition 
region. In the upper inset, we consider the 6 = 10m data 
in the far zone proper. (The other data sets also display 
similar oscillations in their Hamiltonian constraint violations 
in the far zone, though the amplitude of these oscillations 
is too small to be visible on the scale we use to display the 
oscillations of the b — 10m data's Hamiltonian constraint.) 



near zone only accounts for retardation perturbativcly, 
and thus accumulates large phase errors beyond r > A: 
It thus even has a considerable phase error in much of 
our current transition region. As was (briefly) discussed 
in the previous subsection, a resummation of the far zone 
metric might reduce its constraint violations closer to the 
binary, though we did not pursue this here. 

We compared the constraint violations of the new data 
with the old data (for our standard test system) in Fig. Q] 
and Tabic Q] in the introduction, though we deferred a 
more detailed discussion to here: First, it is important 
to realize that the comparison is somewhat misleading, 
since each paper's data are in a different coordinate sys- 
tem. (The data from Paper I arc in the same harmonic 
coordinate system as this paper's data in the near zone. 
However, this is only true perturbativcly in the inner 
zone, where the black hole background is in a coordinate 
system that is not horizon-penetrating.) 

Second, while our data's Hamiltonian constraint vi- 
olations arc not appreciably better than those of the 
data from Paper II, even though we have matched to 
higher order, this is not unexpected: Even though we 
used horizon-penetrating coordinates for the black hole 
metrics, we used standard PN harmonic coordinates for 
the PN metric; these coordinates are singular at the hori- 
zon. While the merged metric has no coordinate singu- 
larities, the PN metric's coordinate singularity increases 



the constraint violations in the transition regions, making 
them comparable to those from Paper II: Paper IPs data 
use a PN metric with no coordinate singularity as well as 
even further resummation of the black hole backgrounds 
than we have employed here, leading to particularly small 
constraint violations. 

Third, the increased momentum constraint violations 
near the hole in the new data, compared with either of 
the old papers' data, come from the x-component: The 
y-component of the new data's momentum constraint vi- 
olations is much smaller than that of either of the two 
previous sets. However, the x-components of their mo- 
mentum constraints vanish. (The cc-component of the 
momentum constraint of Paper IPs data only vanishes 
before the transformation to horizon-penetrating coordi- 
nates; it develops two spikes after that transformation.) 

Finally, if one compares Fig. Q] with, e.g., Figs. 14 
and 17 in Paper II, one notices differences in the behavior 
of the data from Papers I and II close to the hole (and 
inside the horizon). This is because we have generated 
the plot using higher-order finite differencing (fourth or- 
der vs. second order) and a higher resolution (0.002m 
vs. 0.0087«) in computing the constraint violations here 
than we did in first computing them in Paper II. It was 
necessary to do this to accurately resolve the constraint 
violations in the inner zone, since the metric components 
diverge rapidly there. Additionally, we have used the ver- 
sion of Paper IPs metric that is in horizon-penetrating 
coordinates, while Figs. 14 and 17 in Paper II were gen- 
crated using the version of the data without that addi- 
tional coordinate transformation: The transformation to 
horizon-penetrating coordinates introduces further struc- 
ture in the data's constraint violations in both the x- and 
y-components. 
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FIG. 17: The constraint violations of the 0(v 4 ) and 0(i> 5 ) sets 
of data minus those of the full extended set of data. These 
are computed along the :r-axis near hole 1 for our standard 
test binary. 

We can also compare the constraint violations of the 
full extended data with those of the 0{v A ) and 0(v 5 ) 
versions. We do this in the inner-to-near transition in 
Figs. [T7] and [TS] for an equal-mass and 3 : 1 mass ra- 
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FIG. 18: The constraint violations ol the 0(v 4 ) and 0(v 5 ) sets 
of data minus those of the full extended set of data. These 
are computed along the :r-axis in the inner-to-near transition 
region for a binary with b = 10m and a mass ratio of 3 : 1. 
(The more massive hole is on the right.) 
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FIG. 19: The constraint violations of the 0(v 4 ) and full ex- 
tended sets of data. These are computed along the x-axis in 
the near-to-far transition region (and far zone proper) to the 
right of hole 1 for binaries with b = 10m and mass ratios of 
1 : 1 and 3 : 1. (The more massive hole is on the right.) 
In the right-hand panel, we zoom in to show the differences 
in the oscillation of Hamiltonian constraint violations in the 
far zone proper. In the inset, we show the difference between 
the Hamiltonian constraint violations of the 0(v 4 ) and full 
extended versions of the data in the near-to-far transition for 
an equal-mass binary. 
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FIG. 20: The Hamiltonian constraint violations and norm 
of the momentum constraint violations of the versions of the 
inner zone metric with perturbative and full time dependence 
for binaries with 6 = 10m and mass ratios of 1 : 1 and 3:1. 
We display these in the vicinity of both holes, which we shift 
in the 3 : 1 case so their associated point particles lie at 
x = ±5m. [In the figure, x c := [l/(q + 1) - l/2]6.] 




FIG. 21: The differences between the Hamiltonian constraint 
violations and norm of the momentum constraint violations of 
the version of the 0(v 4 ) data without fourth-order octupole 
terms (04_NoQct) in the inner zone and the version including 
those terms (04). These are computed along the x-axis near 
hole 1 for our standard test binary. 



tio binary, respectively, and in the near-to-far transition, 
along with the far zone proper, for both of those bina- 
ries in Fig. 1191 [We consider a unequal-mass binary to 
make the differences between the 0(v 4 ) and 0(v 5 ) ver- 
sions more pronounced: Most of the 0(v 5 ) terms in the 
far zone metric vanish for an equal-mass binary.] 

In the latter plot (Fig. [TT)|) , we do not show the dif- 
ferences between the 0(v A ) and full extended versions' 
Hamiltonian constraint in the far zone proper, as they 
agree up to the level of numerical truncation error. We 
also do not show the differences between the 0(v 5 ) and 
full extended versions of the data. These two sets only 



differ substantially in the far left-hand portion of the 
transition region, and even there the differences are sev- 
eral orders of magnitude less than those between the 
0{v A ) and full extended versions. Additionally, we do not 
show the behavior of the constraint violations to the left 
of the smaller hole (hole 2): It is qualitatively similar to 
their behavior to the right of the larger hole shown here, 
except that for an equal-mass binary, the 0(v 4 ) version 
of the data has smaller Hamiltonian constraint violations 
than the full extended version in that region, and for a 
mass ratio of 3 : 1, there is no oscillation in the tran- 
sition region in the 0(v 4 ) data's momentum constraint 
violations. 
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We compare the constraint violations of the different 
inner zone versions in Fig. [201 for binaries with mass ra- 
tios of 1 : 1 and 3:1. In this plot, we do not include the 
version of the inner zone with no fourth order octupole 
pieces, as the inclusion of those terms does not affect the 
constraint violations in the inner zone proper at a level 
above numerical truncation error. However, these terms 
do have a noticablc effect on the constraint violations in 
the inner-to-near transition, as can be seen in Fig. [2T1 
which compares the 0{v A ) data with and without the 
fourth-order octupole terms in the inner zone. (The fact 
that including the fourth order octupole terms in the in- 
ner zone metric makes a much larger difference in the 
merged metric than in the inner zone metric itself in the 
transition regions suggests that the additional terms that 
have the most significant effect are those in the coordi- 
nate transformation, not those in the tidal fields.) 

It is also interesting to consider how accurate the data 
are for different mass ratios. One finds that the con- 
straint violations do not behave quite as well as might be 
desired in the inner-to-near transition regions as one in- 
creases the mass ratio. This is shown in Figs. [22l and [23l 
which plot the Hamiltonian constraint and norm of the 
momentum constraint for mass ratios of 1:1, 3:1, 5:1, 
and 10 : 1. The worst behavior is that of the momentum 
constraint in the transition region near the more massive 
hole (hole 1), which increases as the mass of that hole in- 
creases. The Hamiltonian constraint also increases with 
mass ratio in the inner zones around both holes. (This 
is due to the inclusion of the full time dependence of the 
tidal fields — the inner zone constraint violations decrease 
with mass ratio if one only uses the version of the inner 
zone metric with perturbative time dependence.) The 
behavior of the other constraint violations is nonmono- 
tonic. When one looks at the near-to-far transition and 
far zone proper one finds much better behavior: The con- 
straint violations decrease with increasing mass ratio in 
all of those regions, except for a slight increase in the 
momentum constraint violations in the far zone proper 
for unequal mass ratios. This is visible for a mass ratio 
of 3 : 1 in Fig. [T9j we do not display the results for higher 
mass ratios, since they are not particularly interesting. 

This behavior in the transition regions is primarily at- 
tributable to the choices we have made for the transition 
functions. For instance, it is possible to choose parame- 
ters so that the momentum constraint violations decrease 
around the more massive hole as its mass increases. This 
can be accomplished by taking Winner, A oc r\, as in Pa- 
pers I and II. However, with this choice, the decrease 
in momentum constraint violations around hole 1 occurs 
at the cost of the aforementioned extreme increase in 
constraint violations around hole 2 as its mass goes to 
zero. It should be possible to combine the two choices for 
Winner, A to obtain better behavior for unequal mass ra- 
tios. However, we have chosen to leave such fine-tuning of 
transition functions to future work, contenting ourselves 
with providing examples of workable transition functions 
here. 
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FIG. 22: The Hamiltonian constraint violations along the x- 
axis for a binary with a separation of 10m and mass ratios of 
1 : 1, 3 : 1, 5 : 1, and 10 : 1 (q := mi/mj). (The more massive 
hole — hole 1 — is on the right, and the less massive hole — hole 
2 — is on the left.) For ease of comparison, we have shifted all 
the data so that the point midway between the two particles 
is at x = 0. [In the figure, x c := [l/(q + 1) - l/2]6.] In the 
two insets, we zoom in to show how the inner zone metric's 
constraint violations vary with q, looking at the region around 
each hole in the inset closest to it. 
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FIG. 23: The norm of the momentum constraint violations 
along the :z;-axis for a binary with a separation of 10m and 
mass ratios of 1 : 1, 3 : 1, 5 : 1, and 10 : 1 (q := mi/rrte). 
(The more massive hole — hole 1 — is on the right, and the less 
massive hole — hole 2 — is on the left.) For ease of comparison, 
we have shifted all the data so that the point midway between 
the two particles is at x = 0. [In the figure, x c := [l/(q + 1) — 
l/2]6.] 



IX. CONCLUSIONS 

A. Summary 

We have constructed approximate initial data for a 
nonspinning black hole binary in a quasicircular orbit. 
This dSata set has uncontrolled remainders of 0(v 5 ) 
throughout the timeslicc (including the far zone), along 
with remainders of 0(v 3 [R/b] 4 , v 5 [R/b] 3 ) in the inner 
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zone. We have verified the scaling of the uncontrolled 
remainders with v by checking that the constraint viola- 
tions decrease at least as rapidly as they should when the 
binary's orbital separation is increased. We constructed 
this data set by asymptotically matching perturbed black 
hole metrics onto a PN metric and creating transition 
functions to smoothly interpolate between the various 
metrics. The resulting data do not assume conformal 
flatness and contain the binary's outgoing radiation, in 
addition to the tidal deformations on the holes. (We 
have included the quadrupole deformations through 1PN 
along with the lowest-order octupole deformations.) 

The results of the matching are given in Sec. IV Gl and 
Eqs. (|B1[) : Sec. IV Gl gives directions for how to put to- 
gether the coordinate transformation necessary to place 
the inner zone metric in the same coordinate system as 
the near zone metric (to the order we have matched). 
Eqs. (|B1[) give explicit expressions for the tidal fields 
we obtained. [We also found that the inner and near 
zone mass parameters — i.e., M and mi — agree through 
at least 0(v 3 ).] These tidal fields are then inserted into 
Detweiler's perturbed black hole metric, given in Cook- 
Schecl coordinates in Eqs. (|3.5p and (|3.6| . to give the 
inner zone metric. The near and far zone metrics are 
given in Eqs. (|4.ip and (|6.6p , respectively. We describe 
our method of computing the binary's past phase evolu- 
tion, needed for the far zone metric, in Sec. IVI Al and 
the specifics of how we put together the various zones' 
metrics in Appendix [D] Workable (though surely not op- 
timal) transition functions that smoothly interpolate be- 
tween the various zones' metrics are given in Sec. IVIII Bl 

We have also constructed an extension of these data 
that is accurate through 0(v 5 ) in the near and far zones. 
In addition, this extension includes other higher-order 
contributions to the temporal components of the near 
and far zone metrics that were readily available in the 
literature. We also calculated the fourth-order octupole 
pieces of the tidal fields (and the associated piece of 
the coordinate transformation) in Appendix IB 21 and in- 
cluded them, as well. (The 1PN correction to the elec- 
tric octupole tidal field is among the terms we calcu- 
lated and add here.) Additionally, we calculated the full 
time dependence of the tidal fields (for times much less 
than the radiation reaction timescale) and included it 
in the inner zone metric. This is discussed in Appen- 
dices [Bj] and |DJ See Appendix [C] for a discussion of how 
we put together the extension to the near and far zone 
metrics, and Eq. (|B3|) for the (polynomial part of the) 
coordinate transformation that accompanies the fourth- 
order octupole pieces of the tidal fields. See Table HT1 (in 
Sec. IVIip for an overview of the different versions of the 
metrics we considered in the paper. N.B.: While we did 
not include this in the table, for the sake of clarity, we 
considered all the versions of the near zone metrics with 
and without background resummation, displaying the re- 
sulting differences in the volume elements [for the 0(v ) 
and full extended versions] in Figs. O and [51 We found 
the full extended (all) data set (including background 



resummation in the near zone) to be the best, overall, 
considering constraint violations as well as the inclusion 
of terms that we expect to improve evolutions. 

In the process of obtaining these data, we have devel- 
oped a method of fixing the matching parameters when 
matching a black hole onto a PN background that dif- 
fers from that presented by Taylor and Poisson [73[ and 
is more automatable. We have also obtained the 1PN 
corrections to the magnetic quadrupole and electric oc- 
tupole for a circular orbit using this method, neither of 
which Taylor and Poisson computed. 

The accurate description of the tidal deformations on 
the holes contained in these data should substantially re- 
duce the high-frequency component of the initial spurious 
radiation; the use of a high-order PN metric should do 
the same for the low-frequency component. (The combi- 
nation of the high-order PN metric, including accurate 
expressions for the trajectories, and the reduced junk 
radiation should also give a much better quasicircular 
orbit — see, e.g., [l7l . fl9j].) If these data do indeed re- 
duce the initial spurious radiation, then they can be used 
to directly quantify the effects of using the conformally 
flat initial data currently employed, as opposed to data 
that include many more of the system's expected physi- 
cal properties. In addition, the waveforms generated us- 
ing these data would be ideal for the construction of hy- 
bri d nu merical relativity/post-Ncwtonian waveforms (as 
in [lO0t | ^ : Since the initial data are directly connected 
to the PN approximation, the PN parameters and phas- 
ing that are input into the initial data should accurately 
describe the subsequent evolution. 

Of course, one needs to evolve the data to see whether 
these putative improvements are indeed realized. We 
have already coded this data set into Maple scripts, 
which were then converted to C code. Both the scripts 
and codes are freely available online at [36| to anyone 
who is interested in evolving or otherwise studying the 
data. Any evolutions of these data will need to use ei- 
ther excision [3(| HH, or the turducken approach [34| : 
The black hole perturbations are not valid all the way to 
the holes' asymptotically flat ends, preventing the use of 
standard puncture methods. Additionally, since the data 
only satisfy the constraint equations approximately, one 
may want to project them onto the constraint hypersur- 
face before evolution. To do th is, one would need to use 
a code such as 45|_ 63|, 1 10 ll Il02| : The standard solver for 
puncture data 103| requires conformal flatness. How- 
ever, evolutions of data that only approximately solve 
the constraints are possible: See, e.g., [57l. Ic35j . 



B. Possibilities for future initial data constructions 
with this method 

With this work, the asymptotic matching method for 
generating initial data for nonspinning binary black holes 
in a quasicircular orbit first essayed by Alvi (4lj and fur- 
ther developed in Papers I and II has been taken to the 
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highest order possible without further development in its 
constituent parts: Including higher multipolc orders in 
the inner zone would require the input of nonlinear black 
hole perturbation theory. Iterating to higher orders in 
v would not only require the higher multipoles, but also 
an explicit expression for the 0(v 6 ) pieces of the purely 
spatial components of the near zone metric. Determin- 
ing the higher-order-in-u pieces of the far zone metric 
would cither require calculation of further contributions 
from the outer integrals in the DIRE approach [6(j, or 
obtaining the far zone metric via matching, following the 
post-Minkowskian approach of [76|. 

Nevertheless, it might still be possible to improve the 
initial data at the (formal) order presented here, as was 
done in Paper II for the data in Paper I. For instance, 
one could contemplate converting the near zone met- 
ric to Cook-Schccl coordinates (or some other horizon- 
penetrating coordinate system) in a neighborhood of each 
black hole. This would regain the complete agreement 
between background coordinates, with no coordinate sin- 
gularity at the horizon in the near zone metric, that was 
found to improve the numerical agreement of the metrics 
in Paper II. One could also further tweak the transition 
functions, though this is not likely to produce any dra- 
matic improvements in constraint violations. However, it 
is possible that different choices for the transition func- 
tions could improve the data's properties in evolutions. 
For instance, in the current near-to-far transition, the 
near zone metric is used (blended with the far zone met- 
ric) for r ^> A, where its perturbative treatment of re- 
tardation leads to large phasing errors. It might thus be 
preferable to transition closer to the binary, even at the 
cost of greater constraint violations. (Resumming the far 
zone metric is a possibility for reducing these constraint 
violations, as was discussed briefly in Sec. I VIII Ah 

The prospects for generalizing this construction to in- 
clude eccentricity or spin are good: The ingredients 
are nearly all readily available in the literature. Such 
generalization — particularly the inclusion of spin — would 
be an obvious next step if this initial data set indeed sig- 
nificantly reduces the spurious radiation. Including ec- 
centricity would be straightforward, though algebraically 
involved, and can be carried out with the ingredients we 
have used here, perhaps supplemented with the results 
from [73| . (The evolution of the binary's phase and sep- 
aration needed in the far zone metric c an b e obtained 
through 3.5PN order using the results of [104| .) 

However, it would only be possible to obtain the gen- 
eralization of these data for a spinning binary to 0(v 2 ), 
formally, while still including all the formal quadrupolc 
pieces in the inner zone. This is true even though the 
generalization of Blanchet, Faye, and Ponsot's me tric to 
include spin (and spin-orbit coupling) is available |l05l ]. 9 



Errata for the potentials and equa tions of motion in [10 
given in footnotes 6 and 10 of Il06ll . respectively. 



as are the expressions for the source multi pole s necessary 
to obtain the matching far zone metric jl07l |. (The bi- 
nary's evolution und er radiation reaction is also known 
through 2.5PN [Tol].) The bottleneck is the available 
tidally perturbed Kerr metric [l08j |. which only includes 
the quadrupolc perturbations. We have seen that knowl- 
edge of the octupole perturbations is necessary to carry 
out the matching of all the formal quadrupole pieces at 
0(v 4 ) [which one would need to do to obtain data that 
are formally valid through 0(v 3 )]. However, one could 
use the same philosophy we did when computing our ex- 
tension and add the higher-order terms in the near and 
far zones without computing the matching inner zone 
terms. 
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APPENDIX A: COMPARING COOK-SCHEEL 
AND PN HARMONIC COORDINATES 

In the quasi-Cartesian form of Cook-Scheel coordi- 
nates [63 , the Schwarzschild metric is 
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[The transformation from Schwarzschild to Cook-Schccl 
coordinates is given in Eq. (|3.4p .j For comparison, the 
Schwarzschild metric in PN harmonic coordinates is 
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where g^ 1 ' = 0. (One obtains PN harmonic coordinates 
from Schwarzschild coordinates by R — V — M, where V 
is the Schwarzschild radial coordinate. This is just the 
spatial part of the Schwarzschild-to-Cook-Scheel trans- 
formation.) Obviously, these coordinates have preserved 
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the coordinate singularity of the Schwarzschild metric in 
Schwarzschild coordinates. 

Note that the purely temporal component of the Cook- 
Scheel version has the same form as in PN harmonic co- 
ordinates, but all the other components are different . In 
particular, the Cook-Scheel version has a nonzero shift 
component as well as a slightly more involved nondiago- 
nal piece of the spatial metric. Explicitly, the differences 
between the PN harmonic and Cook-Scheel versions of 
the Schwarzschild metric components are 

9 oo ~ Soo' = 0; 
cs _ h AM 2 X k 

y ok y k (i? + M)2 R , (A3) 

CS _ har = 16Af 4 X k X t 

(R-M){R + M) 3 R 2 

(Of course, this must be interpreted purely algebraically, 
since we are subtracting components in different coor- 
dinate systems.) The difference of the purely spatial 
components scales as 0([M / R] A ) as M/R — > 0, which is 
clearly higher order in the near zone. However, the dif- 



ference of the spatiotemporal components is 0([M/ R] 2 ), 
which is of the same order as the terms we do keep 
in the spatial metric. And, indeed, we do see the 
first term in this expansion appearing in our coordinate 
transformation — see Sec. IV Gl This shows that we have 
given up on the exact agreement of background metrics — 
which was seen to significantly improve the quality of 
the merged metric in Paper II — in order to have horizon- 
penetrating coordinates. 



APPENDIX B: TIDAL FIELDS 

1. Comparison with Taylor and Poisson's results 

To facilitate the comparison of our expressions for 
the tidal fields with those obtained by Taylor and Pois- 
son [zD], we collect the results of our matching here. 
These include the results of the fourth order octupole 
matching from the next subsection, and are all put to- 
gether using Eqs. (|5.3|) to give explicit expressions for the 
tidal fields about hole 1: 
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Here, for the purposes of comparison with Taylor and 
Poisson's results, we have included the time dependent 
pieces we know in the quadrupole fields (the time de- 
pendence falls into the uncontrolled remainders in the 
octupole fields). These arc ordinarily contained in £ k i 
and B k i, since we usually treat them as independent tidal 
fields. 

Taylor and Poisson give explicit expressions for the 
quadrupole tidal fields for a binary in a circular or- 
bit in their Eqs. (1.10)-(1.14) [and with alternate nota- 
tion in Eqs. (7.25)-(7.29)]. The parts of the quadrupole 
fields that Taylor and Poisson and we have both com- 
puted agree: These consist of the electric quadrupole, 
including its 1PN corrections; the time derivative of the 
electric quadrupole with no corrections; and the mag- 
netic quadrupole with no corrections, all evaluated at 



t = 0. In fact, we can recover all of Taylor and Pois- 
son's expressions for the tidal fields, including the full 
time dependence, if we evaluate our expressions for the 
tidal fields at t = and then make the substitutions 
x k — ► Xk cos u)t+y k smut and y k — > —x k sinujt+y k cosw£. 
(This is, of course, only accurate for times much less than 
the radiation reaction time scale.) Thus, even though we 
are not given the full time dependence directly from the 
matching, we can obtain it from our results, since they 
are true for any point in the orbit. 

While we have computed certain higher-order contri- 
butions to the tidal fields that Taylor and Poisson did 
not, we cannot improve upon the formal accuracy of their 
result for the tidal heating [given through 1PN in their 
Eq. (9.4)]. This is the case even though Poisson gives the 
contribution of the octupole fields to the tidal heating in 
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Eq. (12) of [68]: The 2PN correction to the expression for 
the tidal heating involves the unknown 2PN correction to 
the electric quadrupole (along with the known 1PN cor- 
rection to the magnetic quadrupole and the Newtonian 
piece of the electric octupole). 

We can also check the lowest-order pieces of all the 
tidal fields we found are the expected Newtonian ones: 
The Newtonian pieces of the tidal fields can be com- 
puted independently using Eqs. (5.45), (5.50), and (5.56) 
in [77j , along with the obvious generalization of that ref- 
erence's Eqs. (5.45) and (5.50) for the magnetic octupole, 
viz., = -(3/8)e su {k d lp)s f3 u , where (5 k is given by 

Eq. (5.56b) of [77]]. (To reproduce our results exactly, one 
needs to evaluate all of these expressions for the Newto- 
nian parts of the tidal fields at r\ = 0. Also, since these 
expressions for the Newtonian parts of the tidal fields are 
valid in the rest frame of hole 1, we need to use the rela- 
tive velocity of the holes in calculating (3 k .) Our expres- 
sion for the 1PN correction to the magnetic quadrupole 
can be checked in the extrem e ma ss ratio limit against 
that computed by Poisson in [l09l |: This result is given 
in an appropriate form for comparison in an unnumbered 
equation in Sec. VII E of Taylor and Poisson [73| and 
agrees with our computation. 

2. 1PN corrections to the fourth order octupole 
fields 

Even though we cannot obtain all the inner zone oc- 
tupole pieces at fourth and fifth orders without including 



the hexadecapole fields (since these will give octupolar 
contributions to the nonpolynomial part at these orders), 
it is possible to match only the polynomial parts, and, 
in doing so, read off the 1PN corrections to the octupo- 
lar tidal fields. As an illustration, we shall read off the 
1PN correction to the electric octupole (along with the 
lowest-order piece of the time derivative of the magnetic 
quadrupole) by matching the octupole parts of the poly- 
nomial pieces at fourth order. 



Except for the added algebraic complication of keep- 
ing higher-order multipole terms, this calculation pro- 
ceeds precisely analogously to the fourth order calcula- 
tion involving the quadrupole-and-lower multipoles of the 
polynomial part in Sec. IV Fl The only subtlety that we 
should mention is one that was already present in our 
original fourth order calculation. However, it was not a 
potential source of confusion there because we were com- 
puting both the polynomial and nonpolynomial parts at 
once. Now that we want to compute the polynomial part 
by itself, we need to bear in mind that two nonpolynomial 
pieces can produce a polynomial piece when multiplied 
together (e.g., f is a nonpolynomial piece, but f 2 is a 
polynomial piece). Therefore, contributions to the near 
zone metric involving terms such as n^n^/ S 2 will con- 
tain polynomial pieces, since n\ contains a factor of 1/f, 
and 1/S 2 contributes a factor of r [see Eq. (|4.5[) ] . The 
final results are 
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Continuing on to fifth order to obtain the 1PN correc- cated, but would otherwise proceed as above. At sixth 
tions to the magnetic octupole and time derivative of the order we are not so fortunate. If we tried to carry out 
electric quadrupole would be algebraically more compli- even the quadrupolar part of the sixth order polynomial 
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matching, so as, e.g., to read off the 2PN correction to 
the electric quadrupole, we would be stymied by our lack 
of knowledge of (X^ p ) 43 : The r-times-a-polynomial-in- 

x a pieces of (X^ p ) 4 (which we expect to be present, as 
there have been such terms at all lower multipole orders) 
will contribute to (h^) 6 via the b(l/R) 4 term. 



APPENDIX C: HIGHER-ORDER TERMS 
1. Discussion 

Let us first catalogue all the higher-order near and far 
zone pieces that are readily available in the literature 
(besides those already discussed in Sec. lVI A|) : Blanchet, 
Faye, and Ponsot [58| give an explicit expression for the 
near zone 2.5PN metric with the standard PN order 
counting — i.e., with remainders of 0(v 8 ) in goo, 0(v 7 ) in 
<7ofc; and 0(v e ) in gki- The expression Pati and Will give 
for the far zone metric perturbation in terms of deriva- 
tives of the multipoles [reproduced here in Eq. (16.1a!) ] 
has remainders of at least 0(v 6 ) in all components. Ad- 
ditionally, Pati and Will express the metric in terms of 
the metric perturbation through 3.5PN order (with the 
standard PN order counting) in, e.g., Eqs. (4.2) of [6(ij |. 
Further contributions to the far zone metric perturba- 
tion can be obtained from the full expression for the field 
multipole expansion of the near zone contributions to the 
far zone [Eq. (2.13) in [6(3]], combined with the expres- 
sions for the field multipoles, , in terms of the source 
multipoles, I Q and J®, in Eqs. (4.7) of 0. 

The lowest-order (Newtonian) pieces of all the source 
multipoles are known — see Eq. (Dl) in [6l[. Additionally, 
further PN corrections to the source multipole moments 
than we used in Sec. I VII can be found in various places: 
Higher-order corrections to X can be obtained from the 
system's binding energy, which Blanchet gives through 
3PN in Eq. (170) of [631 ■ Similarly, one can get higher- 
order corrections to J k from the expression Kidder, Will, 
and Wiseman give for the syste m's a ngular momentum 
(through 2PN) in Eq. (2.13b) of [Tlg) . We g ave th e 1PN 
correction to the mass quadrupole in Eq. (|6.2c[) . even 
though we did not need it to construct our 0(v 4 ) data. 
The 1PN correction to the mass octupole can be obtained 
(up to the caveats mentioned below) from the expres- 
sion for the three-index Epstein- Wagoner (EW) moment 
given in Eq. (6.6b) in Will and Wiseman [7l|- This EW 
moment also yields the 1PN correction to the current 
quadrupole. Even further corrections to various source 
multipoles can be obtained from the expressions for the 
two- and four-index EW moments that Will and Wise- 
man give. 

We can use Blanchet, Faye, and Ponsot's results to ob- 
tain initial data in the near zone that are formally valid 
through 0(v 5 ). However, to obtain 0(v 5 ) initial data 
in the far zone, one needs the 0(v 6 ) pieces of the spa- 
tiotemporal components of the metric. As mentioned in 



Scc. lVIII including these is slightly problematic: There is 
an outer integral term in h ok that is 0(v 6 ) with the or- 
der counting used in Sec. lVII and is not known at present. 
This term looks, schematically, like XJ k jr 3 and comes 
from the first term in the general expression Pati and 
Will give for A ok in Eq. (4.4b) of [H. [Here A Q/3 gives 
the contributions of the gravitational field to r Q/3 , the 
"effective" stress-energy pseudotensor. See Eqs. (2.5)- 
(2.7) in [6(| for the explicit definitions.] However, it is 
not present in the explicit expression for A ofe in the far 
zone that Pati and Will give in Eq. (6.5b) of This is 
not a problem, since Pati and Will derive this expression 
specifically for computing the contributions to the near 
zone field from the outer integrals. They have thus used 
the "quick and dirty" rule from their Eq. (4.9) to elimi- 
nate any pieces that do not contribute there to the order 
considered. This piece is one of those eliminated, as it 
will contribute an ^-independent term only through a 
time derivative: 10 Since both X and J k are constants of 
the motion, up to radiative losses, the resulting contribu- 
tion would be of considerably higher order than Pati and 
Will are keeping. However, the "quick and dirty" rule 
is not applicable to the far zone, and we indeed have to 
consider this term. 

While it would be, in principle, reasonably straightfor- 
ward to compute the unknown term, following the proce- 
dure given by Pati and Will in [6(| , the calculation would 
be involved enough that we do not attempt it here. In 
fact, it is possible to argue that we can ignore this term 
entirely: Recalling that, in practice, extra factors of 1/r 
always make a term smaller by at least a factor of v in the 
far zone [see Eq. (|6.3p and the surrounding discussion], 
we can choose to count all powers of 1/r after the first 
as 0(v), disregarding post-Minkowskian considerations. 
This "practical" order counting is justified numerically, 
since the factors of 1/r will make the highest-order near 
zone terms (whose far zone analogues will not be present) 
similarly small when evaluated in the buffer zone. Since 
this is where we stitch the near and far zone metrics to- 
gether numerically, it is thus the only place where we are 
concerned with their agreement, in practice. 

We therefore disregard post-Minkowskian powers of G 
in our counting, giving, e.g., (X/r) 2 = 0{v b ), as opposed 
to 0(v 4 ), as it was before. Nevertheless, this new or- 
der counting still allows us to keep all the outer inte- 
gral contributions that we know (and, indeed, all the 
terms we calculated in Sec. lVip . since we are now keeping 
terms to one order in v higher than before. It also gives 
XJ' k /r 3 = 0(v r ), so we can safely ignore this unknown 
term. Thus, all we need to consistently calculate initial 
data through 0(v 5 ) in the far zone are the 1PN correc- 
tions to the current dipole and the contributions to h ok 



Here 1Z is the artificial radius of separation between the near 
and far zones in the DIRE approach, which, as Pati and Will 
demonstrate in Sec. II I of |6d . cannot appear in any final results. 
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from the mass hexadecapole and current quadrupolc. As 
mentioned above, all of these are easily obtainable. 



Now, even though we can obtain initial data that are 
formally accurate to 0(v 5 ) in the near and far zones (with 
the new far zone order counting) , this docs not utilize the 
0(v e ) and 0(v 7 ) terms in the purely temporal compo- 
nent of the metric that Blanchet, Faye, and Ponsot give 
us. While including these terms does not increase the for- 
mal accuracy of the data, even just in the near and far 
zones, it is nevertheless possible that adding such terms 
will improve the data's quality in practice. Of course, if 
we include terms through 0(v 7 ) in the purely temporal 
component of the near zone metric, it seems desirable 
to also include them in the purely temporal component 
of the far zone metric, and it is possible to do so, up 
to outer integral terms. Here we shall simply neglect 
the outer integral terms we do not know, since there ap- 
pear to be some in h 00 that are 0(v 7 ) even with the 
new "practical" order counting, e.g. , ones that look like 
IX kl h <kt> /r 2 . [See Eq. (6.5a) in [Sfl]. 11 ] Since we are al- 
ready adding pieces without regard to formal accuracy, it 
does not make much sense to go to the trouble of calculat- 
ing these outer integral terms here. This is particularly 
true since we are not even sure that we have all the terms 
we need in the 1PN correction to the mass octupolc. (We 
do know some even higher-order contributions to the far 
zone metric, but do not include those, since we do not 
know the matching near zone contributions.) 



Outline of the calculation 



This uncertainty arises because the EW moments pre- 
sented by Will and Wiseman arc missing any pieces that 
are pure traces in the first two indices: Will and Wise- 
man were only interested in using their results to com- 
pute the gravitational waveform [via their Eq. (2.18)], 
which is transverse and tracefree. For instance, the 
(7 / '4)mim2bS kl term present in the 1PN correction to X kl 
in Eq. (|6.2cp is not present in the expression for I^ w in 
Eq. (4.17) of [ri!| . If not for the missing trace term, these 
two expressions would be identical, up to a surface term, 
as can be seen from their definitions: ip^y and X kl are de- 
fined in Eq. (2.19a) of [zj and Eq. (4.5b) of [U, respec- 
tively. Unlike for X kl , we do not have an independent cal- 
culation of the 1PN corrections to T klp . It is thus possible 
that we are missing terms such as rnlbS kl Xi + (1 <-> 2) 
in which (as demonstrated in the next subsection) 
would contribute to the 1PN correction to T klp and thus 
to the 0(v 7 ) pieces of goo- 



11 This equation only contains the parts of A 00 that contribute 
to the near zone metric. But one can check, starting from the 
general expression for A 00 Pati and Will give in Eq. (4.4a) of |6C| . 
that there are not any lower-order contributions that only appear 
in the far zone. 



For the far zone metric, we first calculate the higher- 
order multipole contributions that need to be added to 
the expression for h a ^ given in Eq. (|6.1ap . The 0(v 6 ) 
and 0(v 7 ) pieces of h 00 are 
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These were obtained by substituting the expressions for 
(in terms of 1® and J®) from Pati and Will's 
Eqs. (4.7) into their Eq. (2.13) (both from [gfj). One 
also needs higher-order contributions to the expression 
for goo m terms of h a/3 . We do not need to add any 
other new terms to the expression for g a /3, since the ex- 
pressions required to obtain the 0(v e ) pieces of c/ofc & n d 
0(v 5 ) pieces of gu are the same as those given previously 
in Eq. (|6.5| . This follows because we already included 
the — (l/2)h 00 term in gofc f° r formal consistency, even 
though it only gives 0(v 6 ) terms in actuality with our 
order counting. The expression we need to obtain goo 
consistently through 0(v 7 ) is 
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taken from Eq. (4.2a) in [60(. 

We also need the 2PN corrections to the mass 
monopole, along with the 1PN corrections to the mass 
quadrupolc, mass octupole, and current dipole. As dis- 
cussed previously, all of these except the 1PN correc- 
tion to the mass octupole are given directly in the liter- 
ature: The corrections to the mass monopole come from 
Blanchet's Eq. (170) 1631 a nd those to the current dipole 
from Eq. (2.13b) in |lldj |: the corrections to the mass 
quadrupole are given in our Eq. ()6.2c|) . However, it is 
possible to obtain a very simple expression for its time 
derivative (up to the caveats mentioned above) in terms 
of the three-index EW moment -f^w [which is given in, 
e.g., Eq. (6.6b) in 0], viz., X klp = . Here the 

equality holds up to surface terms (i.e., ones involving 
TV), which we can neglect here. To obtain this equal- 
ity, take a time derivative of the definition of the mass 
octupole, use the conservation law dpT a P = to write 
doT a0 = —dkT ak , and integrate the result by parts, giv- 
ing an expression that equals 51^^^ up to surface terms. 

We can now antidiffcrcntiate the resulting expression 
for I klp to obtain T klp (up to the caveats mentioned in the 
previous subsection). If we do this for an arbitrary orbit, 
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wc do not need to worry about missing terms due to the 
constant of integration: As we shall demonstrate below, 
any contributions to Ikip that are constant at Newto- 
nian order vanish for a circular orbit. (It also turns out 
that any terms that are time-independent for a circular 
orbit also vanish.) For a generic orbit, it is easiest to 
start from the explicit 2-body reduction Will and Wise- 
man give in Eq. (6.6b). The resulting expression is ex- 
actly what one would expect from the forms of the 1PN 
corrections to the mass monopolc and quadrupole, viz., 
l klp = mizf p (l + v\/2 - m 2 /2b) + (1^2) + 0(v 4 ). 

To prove our claim that any terms in the 1PN correc- 
tion to l klp that are constant at Newtonian order for an 
arbitrary orbit vanish for a circular orbit, we construct 
all such possible terms. To do so, we note that the the 
binary's only (nonzero) vectorial Newtonian constants of 
the motion are its Newtonian angular momentum L and 
Laplacc-Runge-Lenz vector A := p x L — [j, 2 mb. (Here p 
is the momentum of the reduced mass \i := mim2/m.) 
Noting that these are 0(c _1 ) and 0(c~ 2 ), respectively, 
we can write down the only two possible 0(c -2 ) symmet- 
ric 3-index (Cartesian) tensors involving only those two 
vectors (and constants). These will then be the only pos- 
sible 1PN contributions to T klp that are constant at New- 
tonian order. They arc, up to numerical factors (which 
could include contributions of || A||/||Z,|| 2 ), A( k 6 lp ^ and 
\\L\\L^ k 5 lp ^ /m. The first of these vanishes for a circular 
orbit, since A does. The second is odd under time rever- 
sal and is thus inadmissablc, since T klp should be even, 
from its definition, given in Eq. (4.5b) in 60]. 

With all these ingredients, we can put together the far 
zone metric in the same way as we did in Sec. lVIl We thus 
obtain its various components to the same order as we 
are keeping the near zone metric (with the above caveats 
about missing terms), viz., with uncontrolled remainders 
of 0(f 8 ) in the purely temporal component, 0(v 7 ) in 
the spatiotemporal components, and 0(v e ) in the purely 
spatial components. 



APPENDIX D: COMPUTATION OF THE 
METRICS 

Here we detail exactly how the metrics are computed in 
the Maple scripts that were used (along with BAM) to 
compute the constraint violations and create the plots. 
(The scripts themselves, and the resul ting C code that 
Maple outputs are available online at |36|). 



1. Inner zone 

To compute the inner zone metric around hole 1, we 
substitute the tidal fields given in Eq. (|B1[) into the ex- 
pression for Detweiler's perturbed Schwarzschild metric 
in Cook-Schcel coordinates that we give in Eqs. (|3.5p 
and (|3.6p . taking M = mi. We then transform using 



the coordinate transformation given in Sec. IV Gl The 
resulting metric thus includes the 0(v 5 ) terms in the 
purely temporal and spatial components, though these 
do not increase the formal accuracy of the initial data. 
We also do not perform any expansions after substituting 
the tidal fields and performing the coordinate transfor- 
mation, so the final, transformed metric also contains 
various other higher-order-in-w terms. The inner zone 
metric around hole 2 is obtained by the same procedure, 
along with the transformations detailed at the beginning 
of Sec. IVAl 

We have considered three versions of the inner zone 
metric: The first version (contained in D4_NoDct) comes 
directly from the matching performed in Scc.lVland only 
contains the pieces that we were able to match onto the 
near zone metric while including all of the multipolar 
contributions at a given order [i.e., up to octupolar or- 
der through 0(v 3 ) and then only up to quadrupolar or- 
der through 0(v 5 )]. The second version (contained in 
04) also incorporates the results of the fourth-order oc- 
tupole matching carried out in Appendix IB 21 — this in- 
cludes the 1PN correction to the electric octupole, but 
only the polynomial part of the accompanying coordinate 
transformation. The third version (contained in 05 and 
all) adds on the time dependence of all the tidal fields 
(for a circular orbit), obtained in the manner described 
in Appendix IB 21 though it still uses the same coordinate 
transformation as before. (We use the 1PN expression 
for to when substituting for the unit vectors in obtaining 
the full time dependence. We leave off the known higher- 
order corrections to u> here since the expressions for the 
tidal fields we obtained by matching came from using the 
1PN version of co.) 

We calculate the third version by substituting T£ki — > 
(T — t)£ki(0) and similarly for Chip in Eqs. (|3.6|) before 
substituting in the tidal fields (with full time depen- 
dence). These substitutions are necessary because the 
TSki and TCkip terms in Eqs. (|3.6p come from the expan- 
sions of Ski and Chip (with full dependence on V, the ingo- 
ing Eddington-Finkelstein coordinate) about V = 0. We, 
however, are only including the full time dependence on 
t, the near zone time coordinate, in the tidal fields, due to 
our method of obtaining this dependence. These expres- 
sions will thus will only contain the t£ u and tCkip pieces 
(when expanded about t = 0), so wc make the above 
substitutions to retain the linear T-dependence given by 
the matching while not including the linear t-dependence 
twice. (We experimented with including the full time de- 
pendence of the tidal fields using T instead of t and found 
that the constraint violations increased.) 

This method of computing the metric deliberately does 
not include the effects of the full time dependence on the 
spatial variation of the tidal fields' contributions to the 
metric [due to their dependence on V; see the discus- 
sion following Eq. (|3.6|) ]. since these would enter at the 
same order as the unknown time derivatives of the tidal 
fields (second derivatives of the quadrupole fields, and 
first derivatives of the octupole fields). Of course, the 
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terms we are keeping are higher-order as well, but since 
they would enter with explicit factors of t in the multipole 
expansion, they would not be entangled with the explicit 
appearances of unknown time derivatives. (This follows 
because the Schwarzschild metric is time-independent in 
the coordinates we use.) What we have done is equiva- 
lent to repeating the matching we have performed at each 
value of (near zone time) t, up to orbital shrinkage effects, 
which arc higher-order than the terms we are consider- 
ing here. (We also have not attempted to include the full 
time dependence of the coordinate transformation for the 
reasons discussed in Sec. IVIII ) 

2. Near zone 

We compute the near zone metric by substituting the 
trajectories for the point particles [obtained in the man- 
ner discussed in Sec. lVI A] into Blanchet, Faye, and Pon- 
sot's metric [given in Eqs. (7.2) of [58]], including the 
background resummation given in Sec. IVIII Al Here 
there are, again, three versions of the metric, one giv- 
ing 0(v A ) data, one 0(v 5 ) data, and one containing 
the complete 2.5PN metric. Recall that one needs all 
the components through 0(u 4 ) [resp. 0(v 5 )} in addition 
to the 0(v 5 ) [resp. 0(v 6 )} terms in the spatiotemporal 
components in order to obtain 0(v 4 ) [resp. 0(v 5 )} data. 
N.B.: In order to perform background resummation on 
the purely temporal component of the complete 2.5PN 
metric, one needs to also subtract the 0([mi/ri] 3 ) por- 
tion of the expansion of the background, viz., 2mf/rf, 
in Eq. (|8.2[) . All these versions are constructed by trun- 
cating the metric components to the desired order before 
substituting in the trajectories. No expansions are per- 
formed after that substitution, since we do not want to 
drop the higher-order terms that we are keeping in the 
trajectories (discussed in Sec. IVI A[) for conformity with 
the far zone metric. Due to an oversight, we did not 
include the 3PN corrections to the rclativc-to-COM re- 



lation in the 0(i> 4 ) version of the near zone metric. 
3. Far zone 

The calculation of the far zone metric follows the DIRE 
approach, detailed in Sec. IVII differentiating multipole 
moments to obtain the metric. We have the same three 
versions of the far zone metric as for the near zone metric, 
and obtain them in the same manner: We expand to the 
desired order after performing all the substitutions ex- 
cept for (ft, ui, and b (i.e., the contributions that vary due 
to secular radiation reaction effects and are discussed in 
Sec. IVI A|) . We also do not include any of the terms due 
to derivatives acting on b: These would give nonzero con- 
tributions, due to the separation's retarded time depen- 
dence from orbital shrinkage, but these terms are quite 
small, both formally and practically — unlike those due 
to the radiation reaction effects in the phase or undiffer- 
entiated separation — so we neglect them. [For instance, 
the lowest-order contribution due to the nonzero time 
derivative of b in the dki [I kl (u) /r] contribution to 500 is 
formally 0(v 9 ). It is also numerically small: The largest 
contribution has a magnitude of ~ 10 -7 when evaluated 
at the intersection of the x-axis and the inner boundary 
of the near-to-far transition for our equal-mass test bi- 
nary viz., x ~ 20m — see Sec. IVIII B I For comparison, 
the contribution of the uncorrected [m\jr)v\ + (1 <-> 2) 
term in that situation is ~ 10~ 4 .] The 0(v 4 ) version uses 
a slightly different order counting than the 0(v 5 ) and full 
extended versions: As discussed in Sec. IVII we choose to 
keep the outer integral terms — here these are the terms 
that look like (m/r) 2 — in the 0(v 4 ) data due to post- 
Minkowskian considerations, even though those terms 
are 0(v 5 ) if one interprets the Pati-Will order counting 
strictly. However, as mentioned in Appendix [Cj we do 
not know any of the higher-order outer integral terms, 
so we simply drop them in the 0(v 5 ) and full extended 
data. 



[1] F. Pretorius, Phys. Rev. Lett. 95, 121101 (2005), gr- 
qc/0507014. 

[2] M. Campanelli, C. O. Lousto, P. Marronetti, and 
Y. Zlochower, Phys. Rev. Lett. 96, 111101 (2006), gr- 
qc/0511048. 

[3] J. G. Baker, J. Centrella, D.-I. Choi, M. Koppitz, and 
J. van Meter, Phys. Rev. Lett. 96, 111102 (2006), gr- 
qc/0511103. 

[4] F. Pretorius, in Physics of Relativistic Objects in Com- 
pact Binaries: from Birth to Coalescence, edited by 
M. Colpi, P. Casella, V. Gorini, U. Moschella, and 
A. Possenti (Springer, Berlin, 2009), arXiv:0710.1338 
[gr-qc]. 

[5] J. G. Baker, W. D. Boggs, J. Centrella, B. J. Kelly, 
S. T. McWilliams, and J. R. van Meter, Phys. Rev. D 
78, 044046 (2008), arXiv:0805.1428 [gr-qc]. 



[6] J. A. Gonzalez, U. Sperhake, and B. Briigmann, Phys. 
Rev. D 79, 124006 (2009), arXiv:0811.3952 [gr-qc]. 

[7] M. A. Scheel, M. Boyle, T. Chu, L. E. Kidder, K. D. 
Matthews, and H. D. Pfeiffer, Phys. Rev. D 79, 024003 
(2009), arXiv:0810.1767 [gr-qc]. 

[8] S. Husa, J. A. Gonzalez, M. Hannam, B. Briigmann, and 
U. Sperhake, Class. Quantum Grav. 25, 105006 (2008), 
arXiv:0706.0740 [gr-qc]. 

[9] B. Abbott et al. (LIGO Scientific Collaboration), 
arXiv:0711.3041 [gr-qc]. 
[10] L. Lindblom, B. J. Owen, and D. A. Brown, Phys. Rev. 

D 78, 124020 (2008), arXiv:0809.3844 [gr-qc]. 
[11] LISA: Laser Interferometer Space Antenna: 
A Cornerstone Mission for the observation 
of gravitational waves: System and Technol- 
ogy Study Report, ESA-SCI(2000)11, URL 



47 



|http : //www . srl . caltech . e du/lisa/documents/ sts_l . 05 .pcf 



[12] M. Hannam et al, Phys. Rev. D 79, 084025 (2009), 

arXiv:0901.2437 [gr-qc]. [40 
B. Aylott et al, arXiv:0901.4399 [gr-qc]. 

M. Boyle, D. A. Brown, L. E. Kidder, A. H. Mroue, H. P. [41 
Pfeiffer, M. A. Scheel, G. B. Cook, and S. A. Teukolsky, 
Phys. Rev. D 76, 124038 (2007), arXiv:0710.0158 [gr- [42 
qc]. 

L. Lindblom, arXiv:0907.0457 [gr-qc]. [43 
H. P. Pfeiffer, D. A. Brown, L. E. Kidder, L. Lindblom, 
G. Lovelace, and M. A. Scheel, Class. Quantum Grav. 
24, S59 (2007), gr-qc/0702106. [44 
S. Husa, M. Hannam, J. A. Gonzalez, U. Sperhake, 
and B. Briigmann, Phys. Rev. D 77, 044037 (2008), [45 
arXiv:0706.0904 [gr-qc]. 

B. Walther, B. Briigmann, and D. Miiller, [46 
arXiv:0901.0993 [gr-qc]. 

T. Bode, D. Shoemaker, F. Herrmann, and I. Hinder, [47 
Phys. Rev. D 77, 044027 (2008), arXiv:0711.0669 [gr- 
qc]. [48; 

J. A. Gonzalez, U. Sperhake, B. Briigmann, M. Han- 
nam, and S. Husa, Phys. Rev. Lett. 98, 091101 (2007), [49 
gr-qc/0610154. 

G. B. Cook, Living Rev. Relativity 3 ( 2000), URL [50 
|http : //www. livingreviews . org/lrr- 200Q-5| 
E. Gourgoulhon, J. Phys: Conf. Ser. 91, 012001 (2007), 
arXiv:0704.0149 [gr-qc]. [51 
S. Brandt and B. Briigmann, Phys. Rev. Lett. 78, 3606 
(1997), gr-qc/9703066. [52 
P. Marronetti, W. Tichy, B. Briigmann, J. Gonzalez, 
M. Hannam, S. Husa, and U. Sperhake, Class. Quant. 
Grav. 24, S43 (2007), gr-qc/0701123. [53 
P. Marronetti, W. Tichy, B. Briigmann, J. Gonzalez, 
and U. Sperhake, Phys. Rev. D 77, 064010 (2008), [54 
arXiv:0709.2160 [gr-qc] . 

W. Tichy and P. Marronetti, Phys. Rev. D 76, 061502 [55 

(2007) , gr-qc/0703075. 

W. Tichy and P. Marronetti, Phys. Rev. D 78, 081501 [56 

(2008) , arXiv:0807.2985 [gr-qc]. ' 

W. Tichy, B. Briigmann, and P. Laguna, Phys. Rev. D [57 
68, 064008 (2003), gr-qc/0306020. 

W. Tichy and B. Briigmann, Phys. Rev. D 69, 024006 [58 
(2004), gr-qc/0307027. 

M. A. Scheel, H. P. Pfeiffer, L. Lindblom, L. E. Kidder, [59 
O. Rinne, and S. A. Teukolsky, Phys. Rev. D 74, 104006 

(2006) , gr-qc/0607056. [60 

D. Shoemaker, K. Smith, U. Sperhake, P. Laguna, 

E. Schnetter, and D. Fiske, Class. Quantum Grav. 20, [61 
3729 (2003), gr-qc/0301111. 

G. B. Cook and H. P. Pfeiffer, Phys. Rev. D 70, 104016 [62 
(2004), gr-qc/0407078. 

A. Ashtekar and B. Krishnan, Liv- [63 
ing Rev. Relativity 7 (2004), URL 
|http : / /www . l ivingreviews . org/lrr- 2004 -10] 
D. Brown, O. Sarbach, E. Schnetter, M. Tiglio, P. Di- [64 
ener, I. Hawke, and D. Pollney, Phys. Rev. D 76, 081503 

(2007) , arXiv:0707.3101 [gr-qc]. 

N. Yunes, Class. Quantum Grav. 24, 4313 (2007), gr- 

qc/0 611128. 

URL|http: //www .physics . f au. edu/~wolf /Research/harmo 
N. Yunes, W. Tichy, B. J. Owen, and B. Briigmann, 
Phys. Rev. D 74, 104011 (2006), gr-qc/0503011. 
[38] N. Yunes and W. Tichy, Phys. Rev. D 74, 064013 
(2006), gr-qc/0601046. 



[65 



[67 



A. Garat and R. H. Price, Phys. Rev. D 61, 124011 
(2000), gr-qc/0002013. 

J. A. Valiente Kroon, Class. Quantum Grav. 21, 3237 
(2004), gr-qc/0402033. 

K. Alvi, Phys. Rev. D 61, 124013 (2000), gr- 
qc/9912113. 

S. Nissanke, Phys. Rev. D 73, 124002 (2006), gr- 
qc/0509128. 

B. J. Kelly, W. Tichy, M. Campanelli, and B. F. Whit- 
ing, Phys. Rev. D 76, 024008 (2007), arXiv:0704.0628 
[gr-qc]. 

L. Blanchet, Phys. Rev. D 68, 084002 (2003), gr- 
qc/0304080. 

W. Tichy, B. Briigmann, M. Campanelli, and P. Diener, 
Phys. Rev. D 67, 064008 (2003), gr-qc/0207011. 
G. Lovelace, Class. Quantum Grav. 26, 114002 (2009), 
arXiv:0812.3132 [gr-qc]. 

G. Lovelace, R. Owen, H. P. Pfeiffer, and T. Chu, Phys. 
Rev. D 78, 084017 (2008), arXiv:0805.4192 [gr-qc]. 

R. A. Matzner, M. F. Huq, and D. Shoemaker, Phys. 
Rev. D 59, 024015 (1999), gr-qc/9805023. 
P. Marronetti and R. A. Matzner, Phys. Rev. Lett. 85, 
5500 (2000), gr-qc/0009044. 

P. Marronetti, M. Huq, P. Laguna, L. Lehner, R. A. 
Matzner, and D. Shoemaker, Phys. Rev. D 62, 024017 
(2000), gr-qc/0001077. 

U. Sperhake, Phys. Rev. D 76, 104015 (2007), gr- 
qc/0606079. 

M. Hannam, S. Husa, B. Briigmann, J. A. Gonzalez, 
and U. Sperhake, Class. Quantum Grav. 24, S15 (2007), 
gr-qc/0612001. 

S. Dain, Phys. Rev. Lett. 87, 121102 (2001), gr- 
qc/0012023. 

S. Dain, Phys. Rev. D 64, 124002 (2001), gr- 
qc/0103030. 

H. P. Pfeiffer, G. B. Cook, and S. A. Teukolsky, Phys. 
Rev. D 66, 024027 (2002), gr-qc/0203085. 

H.-J. Yo, J. N. Cook, S. L. Shapiro, and T. W. Baum- 
garte, Phys. Rev. D 70, 084033 (2004), gr-qc/0406020. 
B. J. Kelly, W. Tichy, Y. Zlochower, M. Campanelli, 
and B. F. Whiting, in preparation. 

L. Blanchet, G. Faye, and B. Ponsot, Phys. Rev. D 58, 
124002 (1998), gr-qc/9804079. 

S. Detweiler, Class. Quantum Grav. 22, S681 (2005), 
gr-qc/0501004. 

M. E. Pati and C. M. Will, Phys. Rev. D 62, 124015 
(2000), gr-qc/0007087. 

M. E. Pati and C. M. Will, Phys. Rev. D 65, 104008 
(2002), gr-qc/0201001. 

L. Blanchet, Living Rev. Relativity 9 (2006), URL 

http : //www. livingreviews . org/lrr- 2006- 4 

H. P. Pfeiffer, L. E. Kidder, M. A. Scheel, and S. A. 

Teukolsky, Comput. Phys. Commun. 152, 253 (2003), 

gr-qc/0202096. 

M. Hoist, L. Lindblom, R. Owen, H. P. Pfeiffer, M. A. 
Scheel, and L. E. Kidder, Phys. Rev. D 70, 084017 
(2004), gr-qc/0407011. 

T. Bode, P. Laguna, D. M. Shoemaker, I. Hinder, 

F. Herrmann, and B. Vaishnav, arXiv:0902.1127 [gr-qc]. 
Osnsei]ajrdnfitBr/U|pnhtiinll unpublished). 

G. B. Cook and M. A. Scheel, Phys. Rev. D 56, 4775 
(1997). 

E. Poisson, Phys. Rev. Lett. 94, 161103 (2005), gr- 
qc/0501032. 



48 



[69] K. S. Thorne and J. B. Hartle, Phys. Rev. D 31, 1815 [90 
(1985). 

X.-H. Zhang, Phys. Rev. D 34, 991 (1986). [91 
C. M. Will and A. G. Wiseman, Phys. Rev. D 54, 4813 
(1996). 

P. Jaranowski and G. Schafer, Phys. Rev. D 57, 7274 [92 
(1998), gr-qc/9712075. 

S. Taylor and E. Poisson, Phys. Rev. D 78, 084016 [93 
(2008), arXiv:0806.3052 [gr-qc]. 

P. Musgrave, D. P ollney, and K. Lake, GRTensorll soft- [94 
ware (1996), URL http://grtensor.phy.queensu.ca/ 
C. M. Bender and S. A. Orszag, Advanced Mathemati- [95 
cal Methods for Scientists and Engineers (McGraw-Hill, 
New York, 1978). [96 
O. Poujade and L. Blanchet, Phys. Rev. D 65, 124020 
(2002), gr-qc/0112057. 

J. B. Hartle, K. S. Thorne, and R. H. Price, in Black [97 
Holes: The Membrane Paradigm, edited by K. S. 
Thorne, R. H. Price, and D. A. Macdonald (Yale Uni- 
versity Press, New Haven, 1986). [99 
S. Detweiler (private communication). 
R. Rieth, in Mathematics of Gravitation. Part II: Grav- 
itational Wave Detection, edited by A. Krolak (Polish [100 
Academy of Sciences, Institute of Mathematics, War- 
saw, 1997), vol. 41 of Banach Center Publications, pp. [101 
71-74, nearly identical to gr-qc/9603043. 

T. Damour, P. Jaranowski, and G. Schafer, Phys. Rev. [102 
D 62, 084011 (2000), gr-qc/0005034. [103 
J. A. Isenberg, Int. J. Mod. Phys. D 17, 265 (2008), 
gr-qc/0702113. [104 
J. R. Wilson, G. J. Mathews, and P. Marronetti, Phys. 
Rev. D 54, 1317 (1996), gr-qc/9601017. [105 
R. M. Wald, General Relativity (University of Chicago 
Press, Chicago, 1984). [106 
H. Stephani, D. Kramer, M. MacCallum, C. Hoense- 
laers, and E. Herlt, Exact solutions of Einstein's field 
equations (Cambridge University Press, Cambridge, [107 
2003), 2nd ed. 

J. A. Schouten, Ricci- Calculus: An Introduction to Ten- [108 
sor Analysis and its Geometrical Applications (Springer, 
Berlin, 1954), 2nd ed. [109 
L. Blanchet (private communication). 

G. Schafer, Ann. Phys. (N. Y.) 161, 81 (1985). [110 
L. Blanchet, Phys. Rev. D 47, 4392 (1993). 
N. Straumann, General relativity: with applications to 
astrophysics (Springer, Berlin, 2004), rev. ed. 



L. Blanchet and T. Damour, Ann. Inst. Henri Poincare 
Phys. Theor. 50, 377 (1989). 

S. Weinberg, Gravitation and Cosmology: Principles 
and Applications of the General Theory of Relativity 
(Wiley, New York, 1972). 

T. Damour, B. R. Iyer, and B. S. Sathyaprakash, Phys. 
Rev. D 63, 044023 (2001), gr-qc/0010009. 
L. Blanchet and B. R. Iyer, Class. Quantum Grav. 20, 
755 (2003), gr-qc/0209089. 

T. Damour, Int. J. Mod. Phys. A 23, 1130 (2008), 
arXiv:0802.4047 [gr-qc]. 

J. D. Jackson, Classical Electrodynamics (Wiley, New 
York, 1999), 3rd ed. 

L. E. Simone, S. W. Leonard, E. Poisson, and C. M. 
Will, Class. Quantum Grav. 14, 237 (1997), gr- 
qc/9610058. 

J. Pullin (private communication). 

B. Briigmann, W. Tichy, and N. Jansen, Phys. Rev. 
Lett. 92, 211101 (2004), gr-qc/0312112. 

B. Briigmann, J. A. Gonzalez, M. Hannam, S. Husa, 
U. Sperhake, and W. Tichy, Phys. Rev. D 77, 024027 
(2008), gr-qc/0610128. 

P. Ajith et al, Class. Quantum Grav. 24, S689 (2007), 
arXiv:0704.3764 [gr-qc]. 

W. Tichy, A new numerical method to construct binary 
neutron star initial data, in preparation. 

F. Pretorius, in preparation. 

M. Ansorg, B. Briigmann, and W. Tichy, Phys. Rev. D 
70, 064011 (2004), gr-qc/0404056. 

C. Konigsdorffer and A. Gopakumar, Phys. Rev. D 73, 
124012 (2006), gr-qc/0603056. 

H. Tagoshi, A. Ohashi, and B. J. Owen, Phys. Rev. D 
63, 044006 (2001), gr-qc/0010014. 

G. Faye, L. Blanchet, and A. Buonanno, Phys. Rev. D 
74, 104033 (2006), errata in 75, 049903(E) (2007), gr- 
qc/0605139. 

C. M. Will, Phys. Rev. D 71, 084027 (2005), gr- 
qc/0502039. 

N. Yunes and J. A. Gonzalez, Phys. Rev. D 73, 024010 
(2006), gr-qc/0510076. 

E. Poisson, Phys. Rev. D 70, 084044 (2004), gr- 
qc/0407050. 

L. E. Kidder, C. M. Will, and A. G. Wiseman, Phys. 
Rev. D 47, 3281 (1993). 



